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“Modeling  of  Spatial  and  Temporal  Dynamics  in  Biological 

Olfactory  Systems” 

December  2006 

Summary 

The  olfactory  system  is  a  very  efficient  biological  setup  capable  of  odor  information  processing 
with  neural  signals*  The  nature  of  neural  signals  restricts  the  information  representation  to 
multidimensional  temporal  sequences  of  spikes*  The  information  is  contained  in  the  inter -spike 
intervals  in  each  individual  neural  signal  and  in  inter-spike  intervals  between  multiple  signals. 
A  mechanism  of  interactions  between  random  excitaUons  evoked  by  odorants  in  the  olfactory 
receptors  of  the  epithelium  and  deterministic  operation  of  the  olfactory  bulb  is  evaluated  in 
this  project.  Inverse  Frobenius- Perron  models  of  the  bulb's  temporal  sequences  are  fitted  to 
the  interspike  distributions  of  temporally  modulated  receptor  signals.  Ultimately,  such  pattern 
matching  results  In  an  ability  to  recognize  odors,  and  it  offers  a  hypothetical  model  for  signal 
processing  occurring  in  the  primary  stage  of  the  olfactory  system. 
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Chapter  1 

Objectives 


Understanding  the  computational  principles  of  the  biological  neural  circuits  is  a  great  scientific 
challenge  of  the  current  century.  The  olfactory  neural  circuits  are  relatively  well  accessible 
and  studied,  and  are  often  believed  to  be  the  one  that  could  provide  great  insights  and  even 
a  significant  breakthrough  in  our  understanding  of  “How  the  brain  works*. 

Not  only  a  scientific  curiosity  drives  the  Interest  of  studying  the  olfactory  neural  circuits,  but 
also  the  potential  benefits  of  discovering  new  neural  computation  principles.  Indeed,  computa¬ 
tional  abilities  of  the  olfactory  system  are  often  superior  to  the  ones  offered  by  the  engineering 
pattern  recognition  techniques.  Thus,  better  understanding  of  the  olfactory  pattern  processing 
could  open  new  directions  not  only  in  neuroscience,  but  in  computer  science  as  well. 

The  aim  of  the  research  presented  in  this  report  is  to  study  the  principles  of  signal  processing 
occurring  in  the  olfactory  bulb.  A  new  concept  of  the  olfactory  neural  encoding  with  the  inter - 
spike  interval  (ISI)  statistics  has  been  suggested,  and  proposed  to  mimic  the  dynamics  of  the 
olfactory  bulb.  The  Idea  is  that  the  stored  odor  pattern  (ISI  distribution)  corresponds  to  one  of 
the  eigenvectors  of  the  transition  matrix  of  associated  Markov  process. 

A  step  towards  Implementation  of  this  concept  in  a  biologically  realistic  neural  circuit,  a  circuit 
of  integrate-and-fire  neurons,  has  also  been  made.  The  basic  principle  used  in  the  model  is  that 
the  probabilities  of  occurrence  of  particular  Interspike  intervals  in  the  generated  spike  train  are 
controlled  by  a  pattern  of  applied  input  currents.  This  way  the  model  maps  input  patterns  to  the 
ISI  distributions.  It  is  also  suggested  how  such  encoding  with  the  IS!  statistics  could  take  place 
in  biological  olfactory  systems. 

The  Idea  of  controlling  the  dynamics  and  statistics  of  a  neural  dynamical  system  with  an 
external  input  has  been  also  studied  In  a  theoretical  model  of  a  chaotic  system,  based  on  the 
macro-dynamics  of  the  olfactory  bulb  (OB).  The  mechanism  of  the  input-controlled  bifurcation 
in  a  model  of  chaotic  neuron  has  been  proposed.  This  mechanism  enables  the  neural  system  to 
map  the  input  space  to  the  different  modes  of  Its  own  dynamics. 

The  mapping  of  input  patterns  to  the  neural  output  dynamics  has  also  been  studied  in  the 
proposed  model  of  the  olfactory  cortex  dynamics,  where  the  cortical  pattern  recognition  has  been 
investigated.  In  the  developed  model  the  temporal  structure  of  the  olfactory  bulb  input  controls 
spatial  dynamics  of  the  cortical  neural  ensemble. 
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Chapter  2 

Neural  Encoding  in  the  Olfactory 
Bulb 


Biological  olfactory  system  can  be  seen  as  an  engineering  device  working  as  a  pattern  recognition 
machine  with  associative  memory.  It  maps  specific  features  of  the  odor  molecules  detected  by 
the  sensors  of  the  receptor  neurons  in  the  olfactory  epithelium  to  the  inherent  dynamics  of  the 
olfactory  subsystem  called  olfactory  bulb  (OB)  in  mammals  and  antennal  lobe  (AL)  in  insects 
(1.  2.  3.  4.  5.  6.  7].  The  spatio-temporal  patterns  of  this  dynamics  are  then  processed  and 
recognized  by  the  olfactory  cortex  in  mammals  and  mushroom  body  (MB)  in  insects  (8.  9|. 

Variety  of  computational  models  and  abstract  concepts  attempt  to  explain  how  olfactory  sys¬ 
tem  functions.  A  number  of  them  study  the  global  neural  dynamics  of  the  olfactory  bulb,  in 
particular,  the  possible  role  of  chaos  and  bifurcations  [10.  11.  12,  13|  as  well  as  statistics  |14| 
in  the  the  olfactory  information  processing.  More  biologically  oriented  models  explore  the  dy¬ 
namics  of  the  olfactory  bulb  microcircuits  focusing  on  the  coupling  of  neural  oscillators  |15). 
emergence  of  synchronization  [16.  17|  and  the  mechanism  of  pattern  formation  [18.  19|.  The 
recognition  of  these  spatio-temporal  patterns  is  then  studied  in  the  models  of  olfactory  cortex 
[8.  9).  Still,  more  experimental  and  computational  efforts  are  needed  to  clarify  the  basic  infor¬ 
mation  processing  principles  of  the  olfactory  system. 


2. 1  Olfactory  epithelium 

The  receptor  neurons  of  olfactory  epithelium  in  nasal  cavity  transform  certain  properties  of  the 
odor  molecules  to  the  spike  trains  which  they  send  to  olfactory  bulb  (Fig  2.1).  The  odorant 
receptors  on  the  cilia  of  the  sensory  neurons  react  with  the  molecules  having  specific  features.  An 
Individual  sensory  neuron  posses  one  or  a  few  kinds  of  this  receptors,  the  total  number  of  which 
is  about  1000  [20.  2 1 1.  Thus  each  sensory  neurons  is  sensitive  to  one  or  a  few  properties  of  the 
odor  molecule,  although,  the  origin  of  these  molecules’  properties  crucial  for  odor  identification 
is  not  clear  |22[.  The  sensory  neurons  with  identical  receptors  are  spread  uniformly  in  the 
epithelium  in  the  4  zones  [20] .  The  larger  the  concentration  of  the  odor  presented  in  the  nasal 
cavity,  the  more  receptors  are  bound  by  the  molecules,  the  larger  the  number  of  the  sensory 
neurons  are  excited  and  the  larger  are  the  frequencies  of  their  firing.  This  way  epithelium  seem 
to  encode  the  concentration  of  the  odor  molecules  with  the  number  of  activated  sensory  neurons 
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Figure  2.1:  Olfactory  bulb  functional  anatomy.  Abbreviations:  M  -  Mitral  cells:  T  -  tufted  cells; 
Gr  -  Granulle  cells:  PG  -  Periglomerular  cells;  OSN  -  Olfactory  sensory  neurons:(From  Mori  K, 
Nagao  H.  and  Yoshihara.  Y.  Science.  286(5540):71 1-715,  1999). 

and  the  frequency  of  their  firing  (23,  21).  These  spike  trains  from  the  input  which  olfactory  bulb 
receives  from  the  olfactory  epithelium. 


2.2  Functional  anatomy  of  olfactory  bulb 

2.2.1  Glomeruli  layer 

The  axons  (output  cables)  of  the  sensory  neurons  with  the  same  receptor  type  converge  in  the 
OB/AL  to  a  tight  axon  bundle  which  is  called  glomerulus  [20.  6,  21).  Although  the  glomeruli 
might  receive  some  inputs  from  a  different  receptor  types  [20).  it  is  generally  assumed  that 
each  glomerulus  corresponds  to  a  specific  odor  feature  (5).  This  is  the  first  and  basic  signal 
transformation  in  the  olfactory  bulb:  Different  odors  form  specific  spatial  patterns  of  the  excited 
glomeruli,  which  are  activated  by  the  corresponding  receptor  neurons.  Neighboring  glomeruli 
are  interconnected  via  periglomerular  interneurons  (Fig.  2.1).  The  glomeruli  excite  the  adjacent 
periglomerular  cells  and  get  inhibitory  feedback  from  them,  which  modulates  the  glomeruli  ac¬ 
tivity.  It  is  suggested  that  this  modulation  could  be  sharpening  the  spatial  patterns,  which  might 
be  realized  via  wlnner-take-all  competition  [21),  as  in  the  the  visual  system,  although  there  are 
arguments  against  it  (24). 
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2.2.2  Mitral  cells  layer 

Each  glomerulu9  sends  signals  to  one  or  a  few  of  the  specific  mitral  cells  (MC)  In  OB  or  to  the 
projection  neurons  (PN)  in  AL.  MC  and  PN  are  the  principal  (output)  neurons  of  the  OB/AL  (Fig. 
2.1).  The  dynamics  of  mitral  cells  gets  modulated  by  the  interneurons  In  the  similar  way  as  Is 
the  dynamics  of  glomeruli:  Mitral  cells  excite  the  granule  cells  In  return  get  inhibited  by  them 
121). 

Interactions  In  this  excitatory-inhibitory  circuit  supposedly  shape  the  odor-specific  temporal 
patterns  of  the  mitral  cells/projection  neurons  [25.  261.  However,  the  exact  mechanism  of 
formation  of  the  spatio-temporal  patterns  and  their  correlation  with  the  input  odor  patterns  is 
not  clear  [27]. 


2.3  Spatio-temporal  dynamics 

2.3.1  Spatial  patterns 

The  only  clearly  known  encoding  principle  in  the  olfactoiy  bulb  is  that  the  odor  invoked  spatial 
activation  patterns  of  the  glomeruli  and  mitral  cells  of  OB/AL  are  correlated  with  the  odor  com¬ 
ponents  [2,  3,  4,  5,  6|.  Invoked  spatial  patterns  also  proved  to  be  unique  for  different  odors  and 
their  concentrations  (1,  7], 

A  number  of  experiments  show  that  the  neighboring  glomeruli  and  corresponding  mitral  cells 
demonstrate  similar  tuning  properties  [4,  5].  This  suggests  that  the  OB/AL  might  work  as  an 
analog  of  a  Kohonen  map  (28].  There  is  data  that  shows  that  the  mapping  of  the  olfactory  bulb 
may  be  more  complicated  than  that:  Odor  receptor  map  in  the  glomeruli  array  seems  to  have 
hierarchical  and  fractal  dike  structure  II],  It  should  be  noted  that  the  idea  that  spatial  proximity 
of  the  glomeruli  and  mitral  cells  corresponds  to  the  similarity  of  their  tuning  properties  is  not 
universally  accepted  [27],  and  it  is  also  not  clear  how  the  tuning  properties  should  be  defined 
when  a  neuron  participates  in  a  spado- temporal  dynamics  (described  in  the  following  sections). 
Overall,  it  is  widely  believed  that  spatial  acdvlty  patterns  of  mitral  cells  and  glomeruli  are  closely 
related  to  the  odor  encoding,  but  the  exact  structure  of  the  odor -spatial  pattern  correspondence 
is  unclear. 

2.3.2  Spatio-temporal  activity 

Not  only  spatial,  but  also  temporal  patterns  of  activation  proved  to  contribute  to  the  odor  en¬ 
coding  in  the  mammal  olfactory  bulb  |29,  3,  7]  and  insect  antennal  lobe  [30,  31,  26],  As  was 
mendoned  above,  the  neighboring  glomeruli  inhibit  each  other  via  perlglomerular  cells  do  the 
same  with  each  other  via  granule  lnterneurons  (21],  These  interacdons  are  believed  to  be  the 
the  origin  of  the  temporal  dynamics  In  olfactory  bulb.  The  inldal  odor-invoked  spadal  pattern  of 
the  mitral  ceUs  activity  undergoes  certain  transformations  as  some  of  the  mitral  cells,  active  at 
first,  become  inhibited,  while  the  others,  silent/suppressed  at  the  beginning  of  the  odor  onset, 
get  excited  later  (Figs,  2.2  and  2.3}. 

Experimental  data  suggests  two  basic  ideas  about  the  nature  of  spatio-temporal  dynamics. 
Some  of  the  experiments  report  slow  dynamics  of  the  spado- temporal  patterns  In  a  “diffusive" 
manner  |29[,  where  the  Initial  spadal  pattern  does  not  change  significandy:  only  fraction  of 
the  initially  acdve  regions  stops  firing,  and  some  of  the  neighboring  cells  get  recruited  later 
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Figure  2.2:  Activation  patterns  in  the  olfactory  bulb  (From  Spors  H,  and  Grlnvald  A,  Neuron, 
34(2):301-15,  2002.) 


(Fig.  2.2).  It  was  shown  that  such  transformation  optimizes  the  odor  representation  |29). 
The  similar  odors  which  are  close  to  each  other  in  the  odor  space  would  also  produce  similar 
initial  spatial  representations  in  the  mitral  cells  activity.  During  the  transformation,  the  spatial 
activity  patterns  of  the  similar  odors  would  become  more  and  more  distinct,  which,  supposedly, 
facilitates  their  recognition  further  in  the  olfactory  cortex. 

Another,  more  radical  idea  about  the  role  of  spatio-temporal  sequences  in  the  odor  encoding 
is  based  on  experimental  data  on  the  insect  antennal  lobe.  In  the  locust  AL.  an  initial  spatial 
pattern  of  odor  invoked  activity  undergoes  dramatic,  often  cyclic  changes  (30,  26],  in  contrast 
to  comparatively  smaller  changes  in  the  olfactory  bulb  spatial  dynamics  described  above.  This 
data  clearly  shows  that  it  is  not  just  the  initial  spatial  pattern,  but  it  is  the  whole  spatio-temporal 
dynamics  that  encodes  the  odors  [24.  27]  (See  Fig  2.3). 

So,  in  olfactory  bulb  and  antennal  lobe  an  odor  invokes,  at  first,  a  spatial  activation  pattern 
of  the  principal  (mitral  or  projection)  neurons,  which,  then,  transforms  gradually,  in  a  certain 
spatio-temporal  sequence,  which  is  shown  to  be  odor  specific.  However,  what  exactly  the  se¬ 
quence  of  spatial  patterns  and  the  latencies  of  their  appearance  encodes  is  yet  to  be  clarified. 

Some  of  the  hypothesis  suggest  that  the  temporal  structure  of  these  sequences  may  be  re¬ 
lated  to  the  concentration  of  the  odors  [32,  19|.  Indeed,  it  was  shown  that  the  increasing  odor 
concentrations  reduce  response  latencies  of  the  mitral  cells  in  the  odor  specific  sequences,  and 
recruited  new  glomerular  units  [7], 

These  and  other  results  suggest  that  both  odor  identity  and  its  concentrations  are  encoded  by 
both  spatial  and  temporal  structure  of  the  activity  of  mitral  cells.  Interactions  in  the  excitatory- 
inhibitory  mitral-granule  and  glomeruli-periglomeruli  circuit  shape  the  odor-specific  temporal 
patterns  of  the  mitral  cells  [29,  7J.  However,  the  exact  mechanism  of  the  spatio-temporal  pattern 
formation  and  their  correlation  with  the  input  odor  patterns  is  still  unclear  [27]. 
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Figure  2.3:  Odor -specific  spatio-temporal  dynamics  of  the  principal  neurons  in  the  locust  an¬ 
tennal  lobe.  Four  columns  (A-D)correspond  to  the  application  of  four  different  odors.  The  upper 
rows  show  the  spike  train  of  the  projection  neurons  (PN1,  PN2  and  PN3),  and  the  lower  rows 
demonstrate  the  histograms  of  the  spike  trains,  repeatedly  induced  by  the  same  odor  (From 
Laurent  G.  Wehr  M.  and  Davtdowitz  H,  Journal  of  Neuroscience,  16(121:3837-47,  1996.} 

2.4  Encoding  with  transient  synchronization 

In  addition  to  spatio-temporal  dynamics,  another  mechanism  of  encoding,  synchronization, 
proved  to  work  in  olfactory  system.  There  have  been  debates  for  a  long  time  about  the  origin  and 
functional  significance  of  oscillatory  activity  and  its  synchronization  In  the  brain  in  general  and 
in  the  the  olfactory  system  in  particular  [27j.  One  of  the  possible  mechanisms  of  the  emergence 
of  synchronization  is  believed  to  be  based  on  the  interaction  between  excitatory  and  inhibitory 
circuits  (21).  When  excitatory  neurons  (mitral  cells  in  OB  and  projection  neurons  in  AL)  get 
excited  and  start  to  fire  (probably,  incoherently  at  first),  they  activate  common  inhibitory  cells 
(periglomerular  and  granule  cells  in  OB  and  local  interneurons  in  AL).  Inhibitory  cells,  in  their 
turn,  send  their  feedback  simultaneously  to  several  excitatory  neurons  (mitral  cells /projection 
neurons),  and  this  common  inhibition  presumably  synchronizes  these  microcircuits  [33,  21], 

Experimental  data  on  the  insect  olfactory  system  provided  neuroscience  with  one  of  the  clear¬ 
est  demonstrations  of  the  role  of  synchronization  in  neural  encoding  [31,  26],  It  was  shown  that 
information  is  encoded  not  only  by  the  activity  of  the  individual  neurons,  but  also  by  the  fine 
structure  and  mutual  correlations  of  neurons,  in  particular,  by  their  synchronization,  In  the  lo¬ 
cust  antennal  lobe  the  sequences  of  the  transient  synchronizations  of  specific  neural  ensembles 
turned  out  to  be  responsible  for  the  encoding  and  fine  discrimination  (30,  31,  26],  During  the 
stimulus  onset,  specific  groups  of  neurons  synchronize  their  activity  at  specific  cycles  of  their 
activation.  The  timing  of  this  synchronization  and  the  sequence  of  transient  synchronization 
of  different  groups  of  neurons  proved  to  be  odor -specific  and  reproducible  for  a  given  odor  for 
different  trials  [26], 

Moreover,  in  the  experiments  with  the  honeybees  (31)  where  synchronization  was  selectively 
blocked,  but  the  patterns  of  individual  neurons  remained  unchanged,  the  bees  still  could  dls- 
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criminate  between  distinct  odors,  but  were  no  longer  able  to  distinguish  similar  ones.  This 
experimental  data  suggests  that  the  transient  synchronization  phenomenon  in  a  sense,  adds  a 
new  dimension  to  the  input  space.  It  amplifies  the  difference  between  similar  odor  patterns  and 
optimizes  the  stimulus  representation  [27], 


2.5  Modeling  of  the  olfactory  bulb  dynamics 

Two  basic  approaches  could  be  distinguished  in  the  olfactory  bulb  modeling.  One  of  them  deals 
with  the  macrodynamics  of  the  whole  olfactory  bulb,  assuming  that  it  works  as  a  dynamical 
system  [10.  27,  13].  The  microdynamics  approach.  Instead,  studies  the  dynamics  of  the  local 
circuits  and  Interactions  between  mitral  cells.  Interneurons  and  others  OB  components  [15,  18, 
19,  34]. 

2.5.1  Macrodynamics 

Historically,  W.  Freeman  and  co-workers  [10]  were  the  first  to  suggest  that  macrodynamics  of  the 
olfactory  bulb  could  be  seen  in  terms  of  chaotic  dynamics.  Their  experiments  with  the  olfactory 
bulb  of  a  rabbit  support  the  following  concept:  When  there  is  no  odor  applied,  the  macro  state 
of  the  olfactory  bulb  is  wandering  within  a  chaotic  attractor  and  this  may  be  seen  as  the  process 
of  the  solution  search.  Applied  input  {odor)  shifts  the  system  to  the  one  of  its  low-dimensional 
attractors  (“wings")  that  represent  the  recognized  odor  ]10j.  This  concept  inspired  a  number  of 
models  of  chaotic  networks  and  chaotic  neurons  [15.  35.  1 1,  12] 

Recently,  the  idea  that  one  should  analyze  the  olfactory  circuit  activity  in  terms  of  a  dynamical 
systems  was  developed  further.  Based  on  the  experimental  data  on  the  antennal  lobe  dynamics 
]30,  26],  it  was  proposed  that  the  sequences  of  the  spatial  activity  patterns  in  the  locust  antennal 
lobe  can  be  interpreted  and  modeled  as  the  state  of  the  system  visiting  low-dimensional  attractors 
in  the  specific  sequence  and  that  its  trajectory  is  what  encodes  information  [27,  13].  This 
idea  employs  the  concept  of  “winnerless  competition"’:  the  “winning"  attractor  of  the  network's 
dynamics  is  unstable,  so  the  global  state  of  the  system  does  not  get  stuck  there  forever,  but,  after 
being  in  its  vicinity  for  a  while,  moves  on  to  another  attractor  (13],  Such  dynamics  has  been 
also  simulated  in  a  biologically  realistic  neural  circuit  [36]. 

It  should  be  emphasized  that  this  kind  of  dynamics  is  completely  different  from  the  one  of 
the  Hopfield  neural  networks,  where  the  state  is  just  converging  to  the  energy  minimum  from 
the  given  initial  condition.  One  of  the  reasons  behind  it  is  that  In  the  biological  neural  circuits 
input  may  not  be  presented  by  initial  conditions  as  it  is  in  ANN.  The  input  can  be  assumed 
rather  as  a  set  of  the  bifurcation  parameters  that  control  the  system's  dynamics  |15|.  Another 
idea  about  global  dynamics  of  the  olfactory  bulb  was  developed  in  [37].  The  authors  used 
the  Hopfield  neural  network  with  non-symmetric  weight  matrix  to  reproduce  the  experimental 
spatio-temporal  dynamics  of  the  antennal  lobe  [26,  27], 

It  should  be  noted  that  despite  the  efforts  in  modeling  and  theory  of  the  macrodynamics  of 
the  olfactory  bulb,  it  is  still  to  be  shown  that  the  principles  of  dynamical  systems  do  Indeed  apply 
to  the  Information  processing  in  this  system.  More  experimental  data  and  biologically  realistic 
models  are  needed  to  clarify  the  issue  of  the  encoding  with  global  dynamics  in  the  olfactory  bulb. 


1 1 


2.5.2  Dynamics  of  OB  microcircuits 

The  mitral-granule  cell  microcircuit  Is  believed  to  work  as  a  nonlinear  oscillator,  which  is  as¬ 
sumed  to  be  the  origin  of  the  OB  and  olfactory  cortex  oscillations  f  1 5,  21].  One  of  the  first 
models  which  used  coupled  nonlinear  oscillators  to  simulate  olfactory  bulb  dynamics  was  pro¬ 
posed  by  U  and  Hopfield  [34],  Their  model  responded  with  specific  oscillation  modes  to  different 
odors  in  accordance  to  the  Freeman  experimental  data  and  concepts  [10). 

The  modes  of  dynamics  of  a  mitral -granule  pair  of  neurons  in  olfactory  bulb  was  studied  In 
[15].  It  was  shown  that  the  synaptic  strength  of  lateral  (mitral- to-mitral)  connection  may  serve 
as  a  bifurcation  parameter  for  the  mitral-granule  cells  oscillator.  Different  values  of  synaptic 
strength  produced,  via  variety  of  bifurcations,  a  number  of  different  dynamics:  from  fixed  point 
and  limit  cycle  to  strange  attractors  and  chaos. 

Majority  of  the  latest  models  of  olfactory  bulb  focus  on  Its  potential  role  of  segmentation  of 
odor  representation  Into  simpler  patterns  [IS,  19,  38],  This  idea  is  based  on  the  the  concept 
of  temporal  binding  and  segmentation  [39],  which  Is  presumably  at  work  in  the  visual  system. 
The  idea  Is  that  different  components  of  OB  activity  could  be  segregated  In  time  and  processed 
separately. 

The  modeling  efforts  have  been  put  recently  into  studying  the  phenomenon  of  transient  syn¬ 
chronization  in  the  model  of  the  locust  antennal  lobe  [16],  It  was  shown  that  the  mechanism 
behind  this  phenomenon  could  be  based  on  the  Interactions  between  inhibitory  local  neurons 
and  excitatory  projection  neurons,  and  the  competition  between  Inhibitory  circuits.  However,  the 
mechanism  of  the  transient  synchrony  and  its  role  In  the  neural  encoding  is  yet  to  be  clarified. 


2.6  Encoding  with  statistics  of  the  olfactory  bulb  dynamics 

The  characteristic  feature  of  the  information  processing  in  olfactory  system  Is  the  transformation 
of  the  static  odor  signature  (molecular  structure)  Into  the  temporal  sequence  of  neural  activities 
[29,  3,  30,  31,  7,  26],  However,  the  principles  of  such  transformation  are  unknown. 

Most  of  the  computational  models  reproducing  temporal  sequences  of  the  olfactory  bulb/an- 
tennal  lobe  rely  on  the  different  forms  of  neural  competition.  The  main  idea  Is  that  different 
neural  ensembles  win  and  loose  the  mutual  competition  at  different  times  and,  correspondingly, 
get  activated  and  suppressed  in  a  specific  sequence  [36.  18,  19.  34.  38], 

However,  there  Is  an  an  alternative  hypothesis,  which  links  neural  encoding  with  the  statis¬ 
tical  properties  of  the  spike  trains  [40].  In  this  spirit.  It  was  proposed  that  the  odor  features 
(molecules)  could  be  encoded  with  the  probability  density  function  of  the  interspike  Intervals 
produced  by  the  olfactory  bulb  dynamics  [14|.  This  way  the  odors  are  stored  as  the  invariant 
distribution  of  the  probabilities  of  the  interspike  Intervals,  which  corresponds  to  the  one  of  the 
the  eigenvectors  of  the  transition  matrix  of  associated  Markov  process.  So,  the  odors  are  repre¬ 
sented  by  a  dynamical  system  which  produces  corresponding  interspike  interval  distribution. 

This  concept  may  provide  certain  advantages.  The  encoding  with  the  statistics  of  neural 
dynamics  is  invariant  in  time:  the  information  can  be  read  out  from  any  long  enough  part  of  the 
spike  train.  It  may  be  important,  as  the  brains  do  not  really  “know"  when  exactly  the  message  In 
the  spike  train  starts  and  ends. 

The  idea  of  encoding  with  statistics  [14]  is  also  supported  by  the  nature  of  the  biological 
neuron.  The  neuron  is  a  very  unreliable  mechanism:  there  is  no  strict  correlation  between 
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applied  stimulus  and  produced  spike  train.  Thus,  neural  encoding  should  be  based,  at  least  in 
part,  on  the  probabilistic,  rather  than  deterministic  principles  [40). 


2.7  Discussion 

There  has  been  a  significant  progress  in  the  studying  of  the  olfactory  system  recently,  but  still, 
some  major  advances  are  needed  in  the  experimental  measurements  and  the  theoretical  ap¬ 
proaches  to  better  understand  its  functioning.  Experimental  data  provides  some  insights  about 
the  role  of  spatio-temporal  patterns  of  neural  dynamics  in  the  olfactory  bulb  J 1 .  2,  3,  4,  5.  6,  7], 
and  suggests  a  role  of  the  neural  synchrony  in  the  odor  encoding  [30,  27,  33,  31,  26], 

The  first  models  of  the  olfactory  bulb  focused  on  its  global  dynamics,  applying  principles  of 
dynamical  systems  to  the  olfactory  information  processing  [10],  These  ideas  have  been  devel¬ 
oped  further  in  modeling  olfactoiy  [27,  13],  and  other  neural  circuits  jll,  12].  Later,  a  large 
body  of  computational  work  studied  the  formation  of  spatio-temporal  patterns  [36,  18.  19,  38] 
and  different  modes  of  synchrony  [16,  34,  17]  in  olfactory  bulb  microcircuits.  Principles  of  recog¬ 
nition  of  these  patterns  in  the  olfactory  cortex  were  also  explored  [8,  9],  Recently,  a  new  concept 
suggesting  that  the  olfactory  bulb  dynamics  encodes  information  by  means  of  its  statistics  was 
developed  1 14]. 

The  large  number  of  theoretical  approaches  and  computational  models  related  to  olfactory 
system  points  out  that,  probably,  neither  of  them  is  capable  to  explain  fully  the  olfactory  sys¬ 
tem  information  processing.  What  restricts  them  from  being  conclusive  is  the  general  lack  of 
our  understanding  of  some  basic  principles  of  functioning  of  a  single  neuron,  not  to  mention 
the  large  networks.  It  is  not  clear  what  exactly  a  single  neurons  does  in  the  brain,  and  what 
governs  its  dynamics  and  plasticity.  So,  the  future  theoretical  exploration  of  olfactory  system 
is  expected  to  evolve  in  two  directions:  developing  biologically  detailed  models  that  reproduce 
experimentally  measured  dynamics,  and  the  generation  of  new  concepts  about  the  principles  of 
the  global  olfactory  system  dynamics.  The  cooperation  of  these  approaches  should  produce  a 
better  understanding  of  the  olfactory  system  dynamics. 
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Chapter  3 

Signal  Processing  with  Temporal 
Sequences  in  Olfactory  Systems 


Living  organisms  perceive  odors  as  sensations  caused  by  mixtures  of  odorant  molecules.  Such 
molecules  excite  the  olfactoiy  receptors  to  respond  with  increased  activity  which  is  then  passed 
on  to  the  olfactory  bulb  for  detection.  Various  odorant  molecules  excite  different  groups  of 
receptors.  A  superposition  of  these  excitations  constitute  the  odor  as  detected  by  the  olfac¬ 
tory  bulb  [41].  The  relative  concentrations  of  individual  components  constitute  the  odor  type, 
whereas  the  absolute  concentrations  determine  the  odor  intensity.  The  olfactory  bulb  has  the 
task  of  transforming  the  input  obtained  from  the  receptors  into  a  set  of  signals  to  be  interpreted 
by  the  brain. 

The  continuous  quantities,  such  as  molecule  concentrations,  cannot  be  directly  represented 
by  the  signals  produced  by  biological  neurons.  Neurons  produce  spikes  and  only  indirectly  their 
presence  or  absence,  or  time  location  may  carry  continuum  of  information.  The  nature  of  netirai 
signals  is  assumed  to  have  the  following  characteristics: 

1 .  There  is  no  significance  of  the  shape  of  individual  spikes.  They  simply  mark  instances  of 
time  when  the  neurons  fire. 

2.  The  signal  is  a  time  sequence  of  spikes.  Spikes  may  occur  more  or  less  frequently  which 
has  an  effect  on  the  average  value  of  the  signal. 

3.  Spikes  may  occur  in  a  certain  temporal  pattern.  More  precisely,  the  inter- spike  Intervals 
may  follow  a  distinct  and  repetitive  behavior.  This  allows  for  code  division  of  information 
conveyed  by  a  single  signal, 

4.  Two  or  more  signals  may  exhibit  cross-correlation  which  typically  results  from  synchro¬ 
nization  between  the  signal  sources.  If  the  synchronized  signals  assume  a  certain  spatial 
distribuUon,  a  set  of  such  signals  will  manifest  a  spatio-temporal  pattern. 

The  neural  signals  of  the  olfactory  bulb  representing  the  information  about  odors  and  in¬ 
tensities  are  further  interpreted  by  the  brain.  The  olfactory  bulb  functions  as  the  first  signal 
processing  stage.  In  ail  non-blological  designs  the  first  stage  is  responsible  for  the  sensitivity 
and  noise  performance  of  the  entire  detection  system.  The  same  should  hold  true  in  case  of  the 
olfactory  bulb. 
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The  goal  of  this  project  Is  to  identify  the  simplest  method  of  encoding  odor  information  in 
temporal  sequences.  The  input-output  interactions  between  temporal  sequences  can  lead  to  an 
odor  detection  and  encoding  mechanism  in  the  olfactory  bulb. 


3. 1  Temporal  modulation 


The  very  input  of  the  olfactoiy  system,  the  epithelium,  produces  an  enormous  number  of  signals. 
Receptors  are  hard -wired  to  detect  specific  odor  components  and  are  uniformly  distributed  In  the 
epithelium.  The  odor  information  is  therefore  spatially  distributed  across  the  epithelium  and  is 
assumed  to  have  no  temporal  dependency.  Every  odor  and  concentration  can  be  represented 
by  its  “black-and-white  photo"  In  which  the  gray-levels  of  pixels  encode  spiking  activities  of  the 
receptors.  In  this  paper,  the  odor  information  is  assumed  to  be  spatially  distributed  and  static, 
although  there  is  a  strong  evidence  of  various  significant  aspects  of  the  inhale/exhale  rhythm 
and  the  impulse  response  of  the  olfactory  bulb  [42] . 

No  temporal  coding  of  information  is  performed  by  the  individual  receptors.  Simply,  the 
more  molecules  are  present  at  the  docking  sites  of  the  receptors,  the  more  frequent  their  spiking. 
Based  on  response  measurements  and  fitting  of  concentration- response  curves  presented  in  [23J 
and  [431,  the  spiking  frequency  /  of  the  receptors  has  an  asymptotic  dependency  on  the  odorant 
concentration  c  (molarity}.  When  the  odorant  concentration  is  at  the  lowest  detectable  level  c  =  ct, 
the  receptor  fires  at  the  very  slow  rate  of  spontaneous  activity.  When  the  concentration  grows 
infinitely  large,  the  frequency  reaches  saturation  at  the  maximum  firing  rate  of  fm.  The  curve 
/(c)  fits  the  following  definition: 


/  =  — /m  arctan  0(c  -  ct) 


(3.1) 


The  slope  factor  S3  is  expressed  in  terms  of  the  dynamic  range  Ac  defined  as  the  odor  concentra¬ 
tion  at  which  the  frequency  reaches  80%  maximum,  f(ct  +  Ac)  =  0.8/m.  Given  the  dynamic  range, 
the  slope  factor  can  be  determined  as  0  =  tan(0.87r/2)/Ac. 

Concentration  c,  used  in  [231  and  [43],  is  a  logarithmic  quantity  related  to  the  odorant 
molarity  c  =  log  m,  with  m  in  mol/ 1.  The  investigated  odorants  were  anisole,  camphor,  isoamyle 
acetate,  and  limonene,  denoted  by  AN1,  CAM,  ISO,  and  LIM,  respectively.  The  curve  fitting 
resulted  in  the  following  parameters  for  each  odorant  [23]: 


ANI 

CAM 

ISO 

LIM 

/« 

1 1  Hz 

15  Hz 

1 1  Hz 

8  Hz 

Ct 

-6.7 

-8.6 

-7.0 

-7.7 

Ac 

1.1 

1.1 

0.5 

0.3 

Remarkably  linear  curves  are  obtained  if  instead  of  spiking  frequency  /,  the  interspike  in¬ 
tervals  t  =  Iff  are  graphed  versus  the  reciprocal  of  molarity,  referred  to  as  sparsity  s.  Since 
different  odorants  may  have  largely  different  molarity  threshold  levels  ct  =  log mt.  the  reciprocal 
of  the  incremental  molarity  s  =  l/(m  -  mt)  rather  than  the  absolute  value  can  be  used  in  the 
joint  graph  for  various  odorants.  The  parametric  representation  of  relationship  (3.1)  in  the  new 
coordinates  (s,  r)  for  the  introduced  odorants  is  shown  in  Fig.  3. 1 .  The  horizontal  and  vertical 
axes  are  the  incremental  sparsity  s  and  the  interspike  interval  r  expressed  in  terms  of  molarity 
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Figure  3.1:  Interspike  intervals  r  of  receptor  firing  versus  incremental  sparsity  $  of  odorant*  Six 
points  on  each  curve  correspond  to  the  following  logarithmic  concentration  levels  (left-to-right): 
T2Ac,  Ac*  0.8 Ac,  0*6 Ac,  0*4 Ac,  and  0.2Ac  above  ct .  The  axes  units  are  milliseconds  (ms)  and 
megallters  per  moll  (Ml/ mol)* 


m  as  follows: 


s(m) 


1 


m  -  mt 


~fm  arc  tan 
i r 


(3.2) 

(3.3) 


As  can  be  seen  in  the  figure,  diluting  the  odorant  in  the  air  increases  the  Interspike  intervals  at  an 
approximately  constant  rate.  This  may  be  regarded  as  temporal  modulation1  with  the  conversion 
gain  G  —  dr/ds,  which  is  the  slope  of  the  line.  The  left  side  of  each  curve  corresponds  to  the 
receptor  saturation  region.  By  extrapolating  the  curves  to  the  intersections  with  the  vertical 
axis,  a  minimum  interval  r0  for  each  receptor  type  can  be  found  of  value  roughly  around  100  ms. 
This  minimum  interval  may  be  regarded  the  refractory  period  of  the  receptor.  With  just  two 
parameters  r0  and  G  for  each  receptor  type,  the  temporal  modulation  illustrated  in  Fig.  3. 1  can 
be  readily  described  using  first-order  approximation: 


r  =  r0  +  Gs  (3.4) 

The  approximation  can  be  validated  only  within  the  dynamic  range  of  the  receptor,  that  is  outside 
the  saturation  s  >  l/(mt10*e). 


3.2  Odor  characterization  with  interspike  distributions 

An  odor  is  a  superposition  of  a  number  of  basic  odorants.  The  concentration  information  is 
temporally  modulated  at  the  glomerular  Inputs  of  the  olfactory  bulb,  therefore  the  perception  of 

'Term  temporal  modulation  Is  adopted  from  {29] 
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odor  Intensity  must  be  related  to  the  interspike  intervals.  Increasing  the  odor  Intensity  shortens 
the  intervals  at  different  rates  for  each  basic  odorant  due  to  the  differences  in  their  conversion 
gains.  This  provides  some  explanation  why  responses  of  the  mitral  outputs  can  be  different  for 
the  same  odor  at  different  intensities  [441 . 

In  the  glomerular  layer  the  enormous  number  of  inputs  converge  into  much  less  dimensional 
connections  to  the  mitral  cells.  The  glomeruli  are  also  highly  interconnected  between  each  other 
via  periglomerular  interneurons  [45].  Both  inhibitory  and  excitatory  connections  are  present 
within  the  glomeruli  which  indicates  that  a  winner -take-all  mechanism  could  be  involved  before 
the  input  to  the  mitral  cells.  The  presence  of  such  a  mechanism  would  enable  arranging  of  the 
input  interspike  intervals  into  distributions  statistically  representing  the  odor  information. 

Let  R  be  the  number  of  all  types  of  receptors  in  the  epithelium.  This  makes  R  also  the 
number  of  distinct  basic  odorants,  the  basis  for  the  odor  space.  Suppose  the  first  four,  out  of 
all  R  odorants,  are  the  ones  shown  in  Fig.  3.1.  An  odor  at  a  given  intensity  can  be  uniquely 
represented  by  a  vector  s  of  sparsities  of  the  basic  odorants.  For  instance,  an  odor  created 
by  mixing  0.5  mol  of  CAM  and  0.75  mol  of  L1M  with  10OM1  of  air  would  be  represented  by  vector 
s  =  (oo,  50,  oo,  75,  oo, . . . ,  oo)  Ml/  mol.  Vector  2s  would  represent  the  same  mixture  diluted  in  twice 
the  amount  of  air.  In  general,  an  odor,  as  seen  by  the  epithelium,  is  s  =  (si,s2,...,sft).  Terms 
“vector”  and  “basis"  are  understood  to  be  suitable  ways  to  arrange  numbers  rather  than  the 
strictly  defined  terms  used  in  linear  spaces. 

A  much  more  compressed  way  to  describe  odors  is  through  distributions  of  interspike  interval 
probabilities.  This  formalism  may  also  be  more  relevant  to  the  signals  presented  to  the  mitral 
Inputs.  Let  the  interspike  intervals  be  quantized  into  N  ranges  with  cutoff  rmail.  Interval 
is  considered  to  be  a  borderline  between  evoked  and  spontaneous  activity  of  the  receptors.  A 
single  neural  signal  can  represent  an  odor  with  the  interspike  interval  r  probability  distribution 
P-  (pi,pa,  — Pn)-  which  is  formally  a  vector  of  probabilities: 


if  n  <  N 
if  n  =  N 


(3.5) 


The  quantized  representation  of  the  interspike  Interval  distribution  is  chosen  because  it  is  more 
suitable  for  numeric  computations  than  the  probability  density  function.  A  satisfactory  approxi¬ 
mation  of  continuum  can  be  achieved  provided  that  N  is  large  enough. 

Suppose  the  0.5  mol  of  CAM  and  0.75  mol  of  LIM  mixture  with  100  Ml  of  air,  indicated  by  the 
filled  squares  in  Fig.  3.1.  is  presented  to  the  olfactory  epithelium.  Two  kinds  of  receptors  would 
be  activated,  each  responding  with  spikes  separated  by  roughly  90  ms  intervals  and  175  ms  inter¬ 
vals,  respectively,  according  to  Fig.  3.1.  Suppose  also  there  is  twice  as  many  LIM  receptors  as 
those  detecting  CAM.  In  the  winner -take -all  competitions,  the  LIM  receptors  would  have  a  better 
chance  passing  its  signal  compared  to  the  CAM  receptors.  The  described  odor  is  represented  by 
the  filled  bars  of  interspike  interval  probabilities  in  Fig.  3.2.  The  probability  of  the  175  ms  LIM 
intervals  is  twice  the  probability  of  the  90  ms  CAM  intervals. 

Suppose  further  that  the  same  odor  mixture  is  now  diluted  in  twice  the  amount  of  air.  This 
doubles  the  sparsity  of  the  odor,  hence  increasing  the  interspike  intervals  of  both  odorants 
present  in  the  mixture.  The  diluted  odors  are  represented  in  Fig.  3.2  by  the  empty  bars  of  prob¬ 
abilities.  Now  the  LIM  intervals  are  about  220  ms  and  the  CAM  intervals  Increased  to  about  100  ms 
without  a  change  to  the  probability  levels.  Note  that  the  two  odorants  have  different  conversion 
gains  and  modulate  the  temporal  intervals  at  different  rates.  As  the  odor  intensity  changes,  this 
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Figure  3.2:  Odor  composed  of  0.5  mol  of  CAM  odorant  and  0.75  mol  of  LIM  odorant  mixed  with 
100  Ml  of  air  (filled  bars)  and  then  diluted  in  additional  100  Ml  of  air  (empty  bars).  The  bars  have 
width  Tmix/N  and  represent  probabilities  of  respective  time  intervals  as  defined  by  (3.5).  Example 
with  N  =  20  and  rmax  =  350  ms  shown. 

changes  the  probability  pattern.  A  different  hypothesis  of  time  advance  modulation  where  the 
resulting  pattern  is  invariant  under  the  concentration  level  was  introduced  in  (32 j  and  leads  to 
functional  models  [9).  However,  the  neurophysiological  evidence  suggests  that  the  patterns  of 
bulbar  activity  actually  do  change  when  varying  the  odor  intensity  (44). 

The  signal  processing  occurring  between  the  mitral  cells  and  glomerular  layer  is  a  dynamical 
process.  The  information  is  embedded  in  the  time  realizations  of  signals.  It  may  be  retrieved 
only  through  observation  of  these  signals  for  a  period  of  time.  The  probability  distribution  of  the 
interspike  intervals  may  be  retrieved  by  statistically  analyzing  the  neural  signals.  Likewise,  a 
simple  stochastic  process  can  be  modeled  to  have  the  statistical  properties  representing  a  given 
odor  through  the  probability  distribution. 

Suppose,  in  steady-state  after  all  the  transient  response  has  vanished  the  odor  is  represented 
by  the  probability  distribution  p*  defined  according  to  (3.5).  A  Markov  process  (46]  with  the 
Invariant  distribution  equal  p’  could  serve  as  a  first  order  approximation  of  a  dynamical  system 
for  that  odor.  Let  N  x  N  matrix  P  be  the  transition  matrix  of  the  Markov  process 

p(k  +  1)  =  Pp(k)  (3.6) 

Also,  let  the  process  converge  to  p*  in  a  sense  that  p‘  =  lim*-,*,  p{k)  for  almost  every  initial 
distribution  p(0).  The  invariant  distribution  is  the  eigenvector  of  transition  matrix  P  associated 
with  the  unit  eigenvalue:  p*  =  Pp* .  In  this  respect,  the  Markov  process  is  a  dynamical  system  in 
probabilistic  space  S  =  {p  e  [0;  l)w |  52nPn  =  1}  with  a  stable  fixed  point  p\  Further  on,  space  S 
will  be  referred  to  as  the  odor  space. 

Consequently,  an  odor  may  be  associated  with  an  operator  P  :  S  S  in  the  odor  space.  The 
odor  itself  Is  the  stable  fixed  point  of  the  operator.  There  is  a  benefit  of  such  a  representation  of 
odors.  Operator  P  defines  an  odor  indirectly  through  a  definition  of  a  dynamical  system.  It  is  easy 
and  natural  to  generate  realizations  of  neural  signals  using  such  operators,  which  is  suitable  in 
the  modeling  effort.  There  are  many  operators  that  have  the  same  invariant  distribution.  Hence, 
the  same  odor  information  may  be  redundantly  embedded  in  many  different  processes. 

Formally,  a  realization  of  the  introduced  Markov  process  is  a  sequence  of  interspike  Intervals 
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Figure  3.3:  Unit  segment  potential  function.  In  the  total  cost  function,  numbers  Pt e  [0;  1) 
contribute  much  less  than  numbers  P,;  outside  this  range.  Minimization  of  w  will  attract  all  P,/s 
toward  the  Inside  of  the  unit  segment. 


{Tfc}-  Define  the  interval  range  to  be  Tn  —  [^r^;  rmax)  if  n  <  N,  and  =  [rmax;  oc)  otherwise, 
where  the  interval  range  index  n  is  defined  in  the  same  manner  as  in  (3.5).  For  the  sake  of 
modeling,  through  optimization,  a  particular  operator  P  may  be  developed  to  have  p*  as  its 
invariant  distribution  of  interspike  intervals  over  time.  Denote  the  elements  of  the  operator  by 
,  so  that  P  =  [Pjjj.  where  i  and  j  are  the  row  and  column  indices.  Number  Ptj  is  the  probability 
that  in  the  Markov  process  (3.6)  an  interval  from  the  range  T,  will  follow  the  interval  from  the 
range  Tp, 


P  _  PrQfc+1  €  Tj  and  rk  6  Tj) 
lj  Pr (Tk  €  Tj) 


(3.7) 


There  Is  no  closed-form  formula  for  selecting  P^/s  for  a  given  p*.  However,  starting  with  some 
random  P*/ s,  an  optimization  algorithm  can  be  used  to  find  the  Pi/s  as  the  minimum  of  a  suitable 
cost  function.  Since  ail  P</s  are  probabilities ,  they  must  be  numbers  in  the  unit  segment  from  0 
to  1.  This  fact  allows  for  constructing  one  of  the  components  to  be  included  in  the  cost  function, 
namely  the  unit  segment  potential.  For  each  number  Pijf  a  potential  function  v(Pij).  shown 
in  Fig,  3.3  describes  how  distant  Pij  is  from  the  unit  segment: 


vlP  '  _  W)  - 

-  T  .  _  51-6 


(3.8) 


1  +  (2Pij  -  1)- 

Function  v  attains  the  minimum  in  the  middle  of  the  unit  segment  and  is  maximally  fiat  within 
the  segment.  The  maximally  flat  approximation  147)  with  a  rational  function  is  chosen  to  facil¬ 
itate  the  optimization  process.  The  partial  costs  ti(Py)  sum  up  to  the  cost  function  component 
Ev  (P)  responsible  for  keeping  all  the  entries  of  P  within  the  unit  segment: 


N  N 


(3.9) 


i=  1  j=l 


Operator  P  is  a  probabilistic  matrix  in  a  sense  that  all  its  column  vectors  are  normalized 
probability  distributions,  Therefore,  the  column  sums  of  P  must  sum  up  to  1.  Another  cost 
function  component  EC(P)  measures  the  deviation  from  this  requirement: 

n  /  n  \ 2 

Ec(P)  =  J2  U-Epy  (3.i°) 

J=1  \  i=l  / 
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Operator  P  is  a  well  defined  transition  matrix  of  Markov  process  (3.6)  if  EV[P)  +  EC(P)  =  0.  The 
goal  of  the  optimization  procedure  is  to  develop  operator  P  with  the  constraint  that  p*  Is  It's 
principal  eigenvector  associated  with  eigenvalue  1 .  To  simplify  the  operator  synthesis,  matrix  P 
will  be  assumed  to  be  diagonallzable:  P  =  BAB"1.  The  diagonal  matrix  A  =  diag(A)  is  composed 
of  N  eigenvalues  A  =  (At,  A2, . . . ,  Xn)  of  P.  Let  Ai  =  1.  The  convergence  rate  of  the  dynamical  sys¬ 
tem  (3.6)  heavily  depends  on  the  radius  of  the  remaining  A,‘s  for  i  >  1.  Operator  P  is  synthesized 
with  random  Ai's,  for  t  >  1,  with  the  assumption  that  |A<|  <  r  <  1  and  the  radius  r  is  kept  low 
to  improve  the  convergence  rate.  In  the  numerical  experiment  r  was  selected  to  be  equal  to  0.2. 
Operator  P  is  diagonal  in  the  basis  constructed  with  the  column  vectors  of  B.  Since  At  =  1,  the 
first  column  vector  of  B  Is  p\  More  precisely.  B,j  =  p*  for  j  =  1.  All  other  entries  BiJt  for  j  >  1, 
are  variables  in  the  optimization  process.  Their  initial  values  are  selected  randomly  from  the  uni¬ 
form  distribution  in  the  range  (-1;  1).  Pinal  matrix  B  is  found  using  an  optimization  algorithm  to 
minimize  the  cost  function  as  in  the  following  expression: 


min  [EV(P)  +  EC(P)] 


(3.11) 


The  minimized  solution  oftentimes  needs  a  final  touch  to  make  sure  that  P  has  no  negative 
entries,  no  entries  greater  than  one  and  that  column  sums  of  P  are  indeed  1.  This  can  be  done 
by  zeroing  of  negative  values  and  normalization  of  columns.  The  principal  eigenvector  of  P  is 
not  very  sensitive  to  such  trimming  of  P.  Software  package  MATRIX  has  been  developed  (See 
Section  5.1)  in  order  to  generate  a  probabilistic  matrix  according  to  the  principles  described  in 
this  section. 

3.3  Embedding  distributions  in  temporal  sequences 

As  illustrated  in  the  example  shown  in  Fig.  3.2,  the  probabilistic  representation  of  odors  and 
Intensities  fits  well  the  random  nature  of  excitations  received  from  the  olfactory  epithelium.  The 
Markov  model  is  also  a  natural  candidate  for  a  simple  approximation  of  the  dynamics  behind 
spike  interactions  driven  by  the  receptors.  The  olfactory  bulb,  however,  should  be  considered 
a  deterministic  system  which  has  no  random  variables  other  than  the  Input  received  from  the 
epithelium.  Moreover,  the  olfactory  bulb  Is  capable  of  self-excitatory  activity  In  response  to  the 
input.  This  may  be  the  factor  contributing  to  both  high  sensitivity  and  high  selectivity  of  the 
sense  of  smell  [48).  From  this  perspective,  it  seems  reasonable  to  regard  the  olfactory  bulb  as 
an  active  medium  rather  than  a  passive  relay  of  receptor  signals.  The  olfactory  bulb  actively 
produces  firing  activity  in  response  to  the  receptor  signals  [49). 

A  sequence  of  interspike  intervals  complying  with  a  given  interval  distribution  can  be  gener¬ 
ated  in  a  deterministic  dynamical  system.  The  simplest  such  system  is  a  one-dimensional  map 
constructed  by  solving  the  Inverse  Frobenius-Perron  problem  (50).  The  overall  goal  of  the  search 
for  a  sequence  generator  is  to  be  able  to  represent  the  odor  information  by  a  distribution  of  in¬ 
terspike  intervals.  A  simple  shift- map  can  be  constructed  directly  from  the  probabilistic  operator 
used  in  approximation  (3.6),  as  described  in  detail  in  (51). 

First,  a  piece-wise  linear  map  /  :  [0;  Af]  — » (0;  N],  [0;  N]  c  ft  Is  derived  from  probabilities  included 
in  the  operator  P: 
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Figure  3.4:  Example  of  a  piece- wise  linear  shift  map  f(x).  Function  /  is  composed  of  N  contin¬ 
uous  branches  fj  :  [j  —  1;  j)  — *  [0;  TV).  If  x  is  chosen  randomly  from  the  uniform  distribution  over 
the  range  (0;  N),  the  conditional  probability  Ptj  that  /(x)  e  (i  -  1;  i)  given  the  fact  x  e  (j  -  1;  j)  can 
be  evaluated  by  Fy  =  ( | /j- 1  ( [i  -  l;t])||.  Example  with  JV  =  3  shown. 

As  shown  in  Fig.  3.4,  map  /  Is  composed  of  N2  linear  segments  corresponding  to  N2  numbers 
Pij.  The  slopes  of  the  segments  are  simply  1/Fy.  To  evaluate  f(x),  the  pair  of  indices  i  and  j 
appropriate  for  a  given  x  needs  to  be  identified  using  the  condition  provided  by  (3.12). 

By  scaling  the  domain  and  range  of  map  /,  the  dynamical  system  generating  temporal  se¬ 
quence  {xfc }  can  be  defined  with  the  help  of  shift  map  h: 


(3.13) 


Regardless  of  the  initial  condition  chosen,  subsequent  mappings  with  h  will  determine  sequence 
{r*}  whose  distribution  of  values  converges  to  the  invariant  distribution  of  process  (3,6),  A 
deterministic  dynamical  system 


rfc+i  =  h(rk) 


(3.14) 


may  be  regarded  as  a  generator  of  realizations  of  neural  signals  for  a  given  distribution  of  inter  - 
spike  Intervals. 

A  numeric  example  of  shift- map  synthesis  is  shown  in  Fig.  3.5.  Three  interspike  interval 
distributions  p\,  pB,  and  p*c  are  selected  randomly  to  characterize  three  hypothetical  odors  A. 
B,  and  C.  The  bars  represent  probabilities  p„  of  respective  time  intervals  as  defined  by  (3.5). 
The  horizontal  axis  is  normalized  such  that  interval  Nrmax/{N  -  1)  corresponds  to  1.  Shift-maps 
h.A.  hg.  and  he  are  evaluated  for  the  example  odors  and  shown  in  the  middle  row  of  figures 
also  in  time- normalized  coordinates.  The  maps  have  N  disconnected  branches,  however,  vertical 
lines  connecting  the  branches  are  added  to  enhance  the  graphs.  Starting  with  a  random  initial 
interval,  each  map  iterated  3000  times  according  to  (3.14)  produced  a  temporal  sequence.  The 
sequences  are  shown  in  the  bottom  row  of  figures.  Each  interval  in  a  sequence  is  Indicated  as 
a  point  whose  vertical  coordinate  is  the  normalized  time.  This  way  the  density  of  points  reflects 
the  original  distribution  plots  if  rotated  clockwise  by  90°. 
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Figure  3.5:  Synthesis  of  temporal  sequence  generators.  Three  example  interspike  interval  dis¬ 
tributions  with  N  =  20  representing  three  different  odors  are  shown  in  the  top  row.  The  cor¬ 
responding  shift-maps  and  distributions  of  values  of  generated  temporal  sequences  are  shown 
underneath.  The  time  interval  axes  are  normalized  to  the  range  of  (0;  1).  Each  graph  in  the 
bottom  row  contains  3000  points  representing  Interspike  intervals  placed  vertically  according  to 
the  length  of  the  interval. 
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3.4  Frobenius  filter  for  temporal  sequences 

It  Is  broadly  accepted  that  the  olfactory  bulb  provides  support  for  a  pattern  recognition  mech¬ 
anism  for  odor  detection  and  classification.  Not  all  of  the  recognition  is  taking  place  there,  but 
definitely  the  process  is  initiated  in  the  olfactory  bulb.  Assuming  that  the  temporal  sequences  of 
interspike  intervals  are  carriers  for  the  odor  information,  an  implementation  of  signal  processing 
system  (3.14)  can  be  proposed.  Ultimately,  the  goal  is  to  demonstrate  usefulness  of  the  proposed 
mechanism  In  odor  recognition. 

The  proposed  signal  processing  scheme  is  shown  in  Fig.  3.6  and  will  be  referred  to  as  the 
Frobenius  filter.  The  input  to  the  filter  is  a  temporal  sequence  whose  interspike  intervals  are 
determined  by  the  random  variable  r-m  with  values  governed  by  the  probability  distribution  pjn 
defined  as  in  (3.5).  Distribution  pm  characterizes  an  odor. 

The  Frobenius  filter  is  simply  a  shift-map  with  the  feedback  loop  controlled  by  a  random 
switch.  The  switch  operation  is  described  by  a  two-valued  stochastic  process  £  :  {0, 1}  x  N  — *  R. 
The  filter  is  producing  time  intervals  based  on  the  switch  position.  At  every  interval  fc,  the  switch 
position  depends  on  the  value  of  &  governed  by  probabilities: 

Pr{&  =  l)=c  (3.15) 

Pr(&=0)  =  l-e  (3.16) 

where  c  €  [0;  I]  is  a  constant  parameter.  When  this  is  the  case,  the  filter  is  receiving  the  input  rin. 
The  opposite  position  of  the  switch  lets  the  shift-map  determine  the  output  time  interval  based 
on  the  previous  interval  as  in  (3. 14).  The  overall  filter  equation  reads: 

r*+i  =  h[£fcTin  +  {1  ~  (3.17) 

The  notion  of  the  switch  is  an  attempt  to  model  a  competition  between  the  input  and  the  feed¬ 
back.  Its  random  operation  is  inherited  from  the  random  nature  of  the  input  temporal  sequence. 

The  three  shift- maps  Ha,  ha,  and  he,  introduced  in  Fig.  3.5,  are  used  to  Illustrate  the  function 
of  the  Frobenius  filter.  Each  of  the  shift-maps  was  stimulated  at  the  input  by  values  generated 
by  probability  distributions  p\,  p’B,  and  p"c  representing  three  different  odors.  Figure  3.7  shows 
all  possible  input-filter  combinations  arranged  in  the  following  nine  pairs: 

(Pb,/i>i)  (Pp,/i,a) 

(Pi.ftfl)  (Pb^b)  (Pc.  M  (3.18) 

(pj.ta)  (Pfi.ftc)  (Pc.^c) 

In  each  instance,  K  =  20000  values  random  values  were  drawn  from  the  input  distribution  and 
applied  with  probability  c  =  0.5  to  the  filter.  The  values  drawn  were  also  sorted  in  the  ascending 


Pr  =  1— c 

Shift- map 

^  ■) 

Pr  =  c 

Figure  3,6:  Frobenius  filter  is  a  shift-map  with  Input.  Either  the  input  interval  rin  or  the  present 
output  interval  rk  Is  transformed  Into  the  next  output  interval  rk+]. 
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Figure  3.7:  Nine  Instances  of  a  Frobenius  filter  stimulated  by  an  input  distribution  for  20000 
iterations.  Quadratic  distance  di  =  (u,  -  it}2  between  cumulative  distributions  of  interspike 
intervals  at  the  input  u,  and  output  t,  of  the  filter  shown.  The  graphs  are  arranged  according 
to  (3.18). 

order  and  stored.  In  the  sorted  input  sequence  {«,}  the  following  property  holds:  i  <  j  =>  u,  <  Uj. 
When  plotted,  the  graph  of  the  sorted  sequence  would  resemble  the  shape  of  the  cumulative 
distribution  function  of  the  random  variable  r;n. 

The  realization  of  the  sequence  {rfc}  generated  by  the  Frobenius  filter  for  K  iterations  were  also 
sorted  in  the  same  manner.  The  sorted  output  sequence  {t,}  was  then  compared  to  the  sorted 
input  sequence  in  Fig.  3.7.  More  precisely,  the  graphs  in  the  figure  are  the  sequences  of  quadratic 
distances  {{u,  -  i*)2}  in  each  of  the  nine  instances.  The  horizontal  line  is  the  mean  square 
value  of  the  distance  (u*  -  U)2,  As  seen  in  the  figure,  the  input-output  sequences  generated  in 
pairs  (p’B,  hB),  and  (p’c,hc)  are  synchronized  in  a  sense  that  the  quadratic  distance 

between  input  and  output  interval  distributions  is  small.  The  distances  in  all  the  other  pairs  are 
significantly  larger.  By  detecting  low  distance  between  the  input  and  the  output  of  the  filter,  an 
odor  recognition  mechanism  can  be  devised. 

Two  examples  of  realizations  of  the  input  and  output  temporal  sequences  are  shown  in  Fig.  3,8 
and  Fig.  3.9.  The  proposed  mechanism  uses  a  pattern  matching  phenomenon  which  signals 
successful  detection  as  a  decreased  distance  between  parameters  of  the  input  and  the  output 
neural  signals.  The  pattern  matching  is  not  based  on  coherence  of  the  two  signals.  As  shown 
in  the  figures,  no  similarity  in  realizations  of  the  input  and  the  output  can  be  observed  in  either 
the  matched  or  unmatched  odor-filter  pairs.  In  case  of  the  matched  pair,  the  similarity  is  in  the 
statistical  properties  of  the  input  and  the  output  signals. 
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Figure  3.8:  Realizations  of  the  input  (top)  and  output  (bottom)  temporal  sequences  for  a  matched 
pair  (p'A,hA).  A  fragment  containing  100  spikes  shown. 
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Figure  3.9:  Realizations  of  the  input  (top)  and  output  (bottom)  temporal  sequences  for  unmatched 
matched  pair  A  fragment  containing  100  spikes  shown. 

3.5  Discussion 

The  details  of  how  the  cells  of  the  olfactory  bulb  could  encode  the  information  in  the  way  de¬ 
scribed  in  this  project  are  not  discussed  here.  The  goal  set  for  this  work  was  to  describe  the 
simplest  method  of  encoding  information  in  temporal  sequences  and  show  the  input-output  in¬ 
teractions  which  can  lead  to  an  odor  detection  and  encoding  mechanism.  All  the  computations 
are  very  simple.  No  memory  is  necessary,  only  the  last  time  interval  is  locally  kept  in  the  evalua¬ 
tion  of  the  next  time  interval  in  the  output  sequence.  Actual  neurons  are  capable  of  performing 
such  a  storage  with  their  inherent  leaky  integration. 

The  shapes  of  the  shift-maps  shown  in  Fig.  3.5  are  not  crucial  for  in  the  operation  of  the 
Frobenius  filters.  The  proposed  shapes  are  just  samples  of  infinite  possibilities  chosen  for  math¬ 
ematical  simplicity.  In  fact,  any  mapping  that  is  mixing  and  expanding  152]  can  be  used  in 
the  Frobenius  filter.  Such  mappings  develop  continuous  invariant  distributions  and  guarantee 
ergodicity  of  the  temporal  sequence  in  a  sense  that  the  invariant  distribution  is  reachable  from 
any  initial  condition.  Temporal  sequences  at  the  output  of  the  filter  have  a  very  strong  ability  of 
encoding  information  in  the  time  scale,  it  is  sufficient  to  isolate  just  a  few  consecutive  spikes  to 
be  able  to  identify  the  shift-map  which  generated  the  spikes  and  effectively  identify  the  encoded 
odor. 

How  exactly  the  proposed  encoding  concept  may  be  implemented  in  a  biological  neural  circuit 
is  not  directly  addressed  here.  However,  it  would  be  an  intriguing  subject  for  a  possible  future 
research.  A  circuit  of  integrate-and-fire  neurons  can  be  designed  to  produce  spike  trains  with 
controlled  interspike  interval  statistics,  The  software  GENERATOR  developed  for  these  simula¬ 
tions  is  presented  later  in  Section  5.2. 
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Chapter  5 

Appendix 


5. 1  MATRIX  User  Manual 

This  software  generates  a  transition  matrix  of  Markov  process.  The  optimization  technique  pro¬ 
duces  a  corresponding  operator  A  with  a  target  principal  eigenvector  v  associated  with  eigen¬ 
value  1. 

MATRIX  Installation 

1.  Copy  all  files  into  an  arbitrarily  selected  directory  e.g.  c :  /matlab /MATRIX. 

2.  Launch  MATLAB. 

File  list: 


colwise  *  m 
cost .  m 

DiagonalMatrix .m 

eigv,m 

formA .m 

penalty .m 

vector 

Umin 


MATRIX  Tutorial 
Task: 

Use  the  MATRIX  algorithm  to  generate  matrix  A  with  a  target  principal  eigenvector  v  associated 
with  eigenvalue  L 

Solution: 

Step  1,  Choose  the  matrix  size  N: 

»  N-20; 
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Step  2.  Define  global  variables: 

>>  global  v  e 

Step  3.  Store  the  target  eigenvector  v  and  the  array  of  target  eigenvalues  e: 

>>  v=vector (N) 

>>  e=eigv(N) 

Step  4.  Define  initial  matrix  U: 

»  U=rand (N,  N-l) 

Step  5-  Optimize  initial  matrix  U: 

»  Uopt=Umin  <  U  > 

Step  6,  Generate  target  matrix  A: 

»  A=f ormA  C  v , e , Uopt ) 

Step  7;  Save  matrix  A  to  the  file  "matrixA". 

>>  save  matrixA  A  -ASCII 


MATRIX  Function  Reference 

The  following  function  prototypes  are  provided  for  the  advanced  user  These  prototypes  are  useful 
to  the  programmer  if  additional  features  are  to  be  added  to  the  software. 

eolwise.m 
USAGE:  colwise(m) 

Evaluates  the  quadratic  distance  between  1  and  the  column  sum,  for  each  of  the  columns  m  of 
the  matrix  A. 

cost .m 

USAGE:  cost  (U) 

Matrix  cost  function,  which  defines  the  optimization  process  of  the  matrix  A. 

DiagonalMat rix .m 
USAGE:  Diagonal-Matrix  (e) 

Creates  a  diagonal  matrix  with  specified  column  of  eigenvalues  e  on  its  diagonal, 

eigv . m 

USAGE:  eigv  (N) 

Creates  a  column  of  eigenvalues  e  of  dimension  N,  where  the  first  one  is  equal  to  1  and  the  rest 
are  equal  to  0.2, 

formA.m 

USAGE:  formA(v,e,U) 

Creates  matrix  A  with  specified  eigenvectors  v  and  corresponding  eigenvalues  e. 
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penalty .m 
USAGE:  penalty  (x) 

Penalty  function  imposing  constraints  on  the  matrix  A. 

vector ,m 
USAGE:  vector  (N} 

Creates  a  random  normalized  vector  of  dimension  N 
Umin . m 

USAGE:  Umin  (U) 

Optimizes  the  Initial  matrix  U  in  order  to  create  matrix  A  with  minimal  cost  function. 


The  code  of  the  corresponding  files  can  be  downloaded  from: 

http :  /  /ci.  uofl.edu  /cu  rrentwork/lozowskl/  depscorO  1  /report .  final  /MATRIX/ 


5.2  GENERATOR  User  Manual 

GENERATOR  Tutorial 

This  software  generates  a  transition  matrix  of  Markov  process.  The  optimization  technique  pro* 
duces  a  corresponding  operator  A  with  a  target  principal  eigenvector  v  associated  with  eigen¬ 
value  1. 

GENERATOR  Installation 

1.  Copy  all  files  Into  an  arbitrarily  selected  directory  e.g.  C:  /mat lab /generator. 

2.  Launch  MATLAB. 

File  list: 


diagr .m 
f  ilterl .m 
isi  .m 
prob .m 


GENERATOR  Tutorial 
Task: 

Use  the  GENERATOR  algorithm  to  produce  a  spike  train  with  the  input-controlled  interspike 
Interval  distribution. 

Solution: 

Step  1.  Set  an  example  input  pattern:  Input5=[0.41,  0.28,  0.41,  0.31,  0.34): 
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»  Input 5= [0 -41;  0.28;  0.41;  0*31;  0.34]; 

Step  2*  Evaluate  and  plot  the  normalized  input  pattern  InputSnorm: 

>>Input5norm  =  Input5/sum { Inputs) 
ans  = 

0.2343 

0.1600 

0.2343 

0,1771 

0.1943 


>>  bar ( InputSnorm) 


T  !  I  «  t 


Step  3,  Generate  the  corresponding  spike  train  train  of  the  length  T  seconds 

»  T=1Q ; 

>>  train=prob (T, Inputs J ; 

Step  4.  Plot  the  diagram  of  the  interspike  interval  distribution  of  the  generated  spike  train: 

»  bar (diagr {train f Input 5) } 
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GENERATOR  Function  Reference 

The  following  function  prototypes  are  provided  for  the  advanced  user.  These  prototypes  are  useful 
to  the  programmer  if  additional  features  are  to  be  added  to  the  software. 

p  rob .  m 

USAGE:  prob (T, Input5) 

Creates  a  spike  train  of  length  T,  controlled  by  the  input  vector  Inputs.  In  the  output  array  the 
“1"  and  “0"  represent,  respectively,  the  presence  or  absence  of  a  spike  in  a  1ms  interval. 

isi  ,m 

USAGE:  isi  (train) 

Produces  an  array  of  the  interspike  Intervals  of  the  spike  train  train  created  by  the  prob.m  func¬ 
tion. 

diagr  .m 

USAGE:  diagr  (train.  Input) 

Creates  the  diagram  of  the  probabilities  of  occurrence  of  the  the  corresponding  Interspike  Inter¬ 
vals  in  the  generated  spike  train  train  controlled  by  the  input  pattern  Input. 

f  i  Iter 1 .ra 
USAGE:  filterl  (u) 

Realizes  the  winner-take-all  competition:  limits  the  number  of  firing  neurons  in  the  neuron  array 
u  during  one  cycle  to  one. 

The  code  of  the  corresponding  files  can  be  downloaded  from: 
http://ci.uofi.edu/currentwork/lozowski/depscor01/report.final/GENERATOR/ 
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Abstract.  This  paper  proposes  temporal-to-spatia!  dy¬ 
namic  mapping  inspired  by  neural  dynamics  of  the 
olfactory  cortex.  In  our  mode!  the  temporal  structure  of 
olfactory-bulb  patterns  is  mapped  to  the  spatial  dynam¬ 
ics  of  the  ensemble  of  cortical  neurons.  This  mapping  is 
based  on  the  following  biological  mechanism:  while 
anterior  part  of  piriform  cortex  can  be  excited  by  the 
afferent  input  alone,  the  posterior  areas  require  both 
afferent  and  association  signals,  which  are  temporally 
correlated  in  a  specific  way.  One  of  the  functional  types 
of  the  neurons  in  our  model  corresponds  to  the  cortical 
spatial  dynamics  and  encodes  odor  components,  and 
another  represents  temporal  activity  of  association-fiber 
signals,  which,  we  suggest,  may  be  relevant  to  the 
encoding  of  odor  concentrations.  The  temporal-to- 
spatial  mapping  and  distributed  representation  of  the 
model  enable  simultaneous  rough  cluster  classification 
and  fine  recognition  of  patterns  within  a  cluster  as  parts 
of  the  same  dynamic  process.  The  model  is  able  to 
extract  and  segment  the  components  of  complex  odor 
patterns  which  are  spatiotemporal  sequences  of  neural 
activity, 


1  Introduction 

Flexible  object  recognition,  feature  binding  and  segmen¬ 
tation,  attention  focusing,  and  other  pattern  processing 
tasks  are  hardly  handled  by  computational  techniques 
based  on  stationary  principles.  On  the  other  hand,  they 
are  successfully  resolved  by  biological  neural  systems, 
where  not  only  spatial,  but  also  different  kinds  of 
temporal  dynamics  and  correlation  are  believed  to  be  the 
underlying  principles  of  .the  brain’s  abilities  (Fujii  et  at. 
1996;  Malsburg  1992). 

Olfaction  is  an  example  of  such  a  system  in  which 
spatiotemporal  dynamics  has  been  the  subject  of  both 
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numerous  experimental  studies  (Duchamp- Viret  et  al. 
1996;  Joerges  et  al.  1997;  Ressler  et  al.  1994;  Scarda  and 
Freeman  1987;  Wehr  and  Laurent  1996)  and  theoretical 
modeling  (Haberly  and  Bower  1989;  Hendin  et  al.  1998; 
Hoshino  et  al.  1998;  Wilson  and  Bower  1992).  However’ 
most  odor  recognition  techniques  do  not  make  use  of 
temporal  encoding  and  processing.  In  static  systems 
patterns  are  recognized  by  the  stationary  pattern  rec¬ 
ognition  methods,  which  do  not  appear  to  be  linked  to 
biological  temporal  dynamics.  The  lack  of  such  dy¬ 
namics  and  dynamic  processing  principles  in  pattern 
recognition  systems  is  one  of  the  reasons  for  their  poor 
computational  abilities  in  comparison  to  biological 
neural  circuits. 

On  the  other  hand,  there  are  a  number  of  olfactory 
cortex  models  which  are  based  strictly  on  the  system’s 
neuro  bio  logical  description  and  produce  spatiotemporal 
dynamics  similar  to  the  experimental  data  (Ballain  et  al. 
1998;  Wilson  and  Bower  1992).  Although  such  models 
do  not  clarify  the  functional  significance  of  this  dy¬ 
namics,  they  provide  insights  about  its  possible  role  in 
cortical  information  processing. 

Indeed,  biological  data  demonstrate  that  pyramidal 
cells  are  principal  neurons  of  olfactory  cortex  that  re¬ 
ceive  and  integrate  four  major  types  of  input  signals: 
from  the  olfactory  bulb  (OB)  through  afferent  fibers, 
and  also  from  three  functional  cortical  areas  via  asso¬ 
ciation  fibers.  These  four  signals  are  delayed  differently 
by  different  branches  of  fibers,  and  arrive  lo  the  area  of 
convergence  at  different  times,  Moreover,  their  signals 
need  different  times  to  propagate  through  the  dendrites 
and  reach  the  cell  body.  Temporal  correlation  of  in¬ 
coming  afferent  and  associative  signals  proved  to  be 
crucial  for  the  signal  integration  by  a  pyramidal  cell  (see 
Sect.  6.1  for  details)  (Haberly  1998).  We  explore  this 
type  of  temporal  correlation  in  our  model  and  suggest 
that  it  could  give  some  cues  to  understand  how  a  rec¬ 
ognition  of  multicomponent  odors  and  their  mixtures  is 
realized  in  the  olfactory  cortex. 

In  the  ensemble  of  spiking  neurons  described  in 
Sect.  3,  the  neurons  at  different  layers  have  different 
functional  properties.  The  dynamics  of  the  neurons 
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which  are  sensitive  to  the  odor  components  represents 
the  activity  of  the  neural  pool  corresponding  to  the 
spatial  pattern  of  activity  induced  by  this  component 
(Haberly  1998),  Although  another  type  of  neuron,  sen¬ 
sitive  to  the  concentrations  of  the  odor  components,  also 
proved  to  have  its  biological  analog  in  the  frog  olfactory 
cortex  (Duchamp- Viret  et  al.  1996),  its  function  can  be 
interpreted  more  generally.  The  activity  of  this  neuron 
delivers  OB  signals  to  the  cortical  neurons  with  different 
delays,  so  it  may  correspond  to  the  dynamics  of  asso¬ 
ciation  fibers. 

This  ensemble**  architecture,  due  to  its  distributed 
and  dynamical  representation  of  memory,  provides  the 
basis  for  solution,  of  the  coarseness-sensitivity  flexibility 
problem.  A  coarse-enough  system  cannot  distinguish 
fine  variations  of  the  patterns  within  a  cluster.  On  the 
other  hand,  a  sensitive-enough  system  is  not  able  to 
detect  what  cluster  these  slightly  different  patterns  be¬ 
long  to.  The  tempo ral-to-spaiial  mapping  and  distrib¬ 
uted  representation  of  the  model  enable  simultaneous 
rough  duster  classification  and  fine  recognition  of  pat¬ 
terns  within  a  duster  as  parts  of  the  same  dynamic 
process. 

Another  group  of  tasks  handled  by  biological  systems 
includes  feature  binding,  segmentation,  attention  focus¬ 
ing,  and  other  multi  pattern  recognition  problems.  There 
are  experimental  data  that  suggest  that  in  the  brain  these 
tasks  are  solved  with  temporal  processing  (Fujii  et  al. 
1996;  Jinks  and  Laing  1999),  and  there  are  models  that 
propose  possible  mechanisms  (Ambros-Ingerson  et  ah 
1990;  Campbell  and  Wang  1998;  Grossberg  1999;  Hen- 
din  et  al.  1998;  Lysetskiy  et  al.  2001;  Malsburg  1992).  In 
Sect.  4  we  show  that  the  temporal  structure  of  our  model 
allows  us  to  realize  this  muitipattern  processing. 

In  Sect.  2  we  introduce  the  biological  olfactory  sys¬ 
tem  and  its  functional  properties.  Sect.  3  describes  the 
basic  block  layout  of  the  model  and  its  spatiotemporal 
dynamics  which  ensures  the  recognition  flexibility.  An 
example  of  temporal  segmentation  of  odor  patterns  in  a 
mixture  is  presented  in  Sect.  4.  Methods  and  parameters 
of  the  simulation  are  shown  in  Sect.  5,  which  is  followed 
by  discussion  and  conclusions  in  Sect.  6. 


2  Olfactory  system 
2.1  Olfactory  bulb 

An  odor  identity  is  defined  by  a  group  of  physical  and 
chemical  parameters  of  the  odor’s  constituent  molecules 
and  their  relative  concentrations.  However,  these  pa¬ 
rameters  are  not  clearly  determined,  nor  is  the  correla¬ 
tion  between  their  candidates  and  the  odor  properties 
(Wise  et  al.  2000).  For  the  sake  of  simplicity  we  assume 
that  one  constituent  molecule  possesses  one  of  these 
crucial  parameters  and  corresponds  to  one  of  the  odor 
components.  The  odors  can  therefore  be  presented  by 

the  concentration  vector  C  =  _ tcn },  where  cj  is 

the  concentration  of  yth  molecule. 

In  the  olfactory  system,  odors  are  first  perceived  in  the 
olfactory  epithelium  by  different  types  of  receptor 


neurons.  These  neurons  are  sensitive  to  differem  kinds  of 
molecules  and  respond  selectively  to  their  presence  with 
oscillatory  firing.  An  odor  is  further  encoded  into  a 
quasiperiodic  temporal  sequence  of  spatial  patterns  of 
synchronized  oscillatory  neural  activity  of  the  olfactory 
bulb  or  antennal  lobe  (Hoshino  et  al.  1998;  Laurent  1996; 
Scarda  and  Freeman  1987;  Wehr  and  Laurent  1996). 

The  spatial  patterns  of  the  olfactory  bulb/antennal 
Jobe  have  been  proved  to  be  correlated  with  the  odor 
components  (Joerges  et  al.  1997;  Laurent  1996;  Ressler 
et  al.  1994),  but  the  functional  significance  of  their 
temporal  structure  is  unclear.  There  are  two  major  hy¬ 
potheses  of  its  possible  role.  Experiments  such  as  these 
of  Laurent  et  al.  (1996)  and  Wehr  and  Laurent  (1996) 
show  that  the  temporal  structure  of  firing  of  different 
ensembles  contributes  to  encoding  of  odor  identity  in  a 
certain  combinatorial  way.  Another  idea  is  that  tempo¬ 
ral  dynamics  of  the  olfactory  bulb  encodes  the  concen¬ 
trations  of  the  odor  components.  According  to  the 
concept  that  the  odor  components  are  encoded  as  dy¬ 
namic  attractors  of  neural  activity  (Scarda  and  Freeman 
1987),  the  order  in  which  the  state  of  the  system  visits 
these  attractors  and  the  time  the  state  spends  wandering 
around  them  could  be  correlated  with  the  concentrations 
of  the  odor  components  (Hoshino  et  al.  1998).  There  is 
also  experimental  and  theoretical  support  for  the  idea 
that  a  component’s  concentrations  can  be  encoded  by 
precise  timing  of  a  spike,  a  burst,  or  the  phase  of  the 
periodic  firing  of  corresponding  neuron  or  ensemble 
(Duchamp- Vi  ret  et  al.  1996;  Hopfield  1995).  These  hy¬ 
pothesis  do  not  necessarily  contradict  each  other;  they 
could  coexist  and  complement  one  another  or  be  the 
parts  of  a  more  complicated  neural  coding  scheme. 

Our  model  realizes  a  mapping  of  the  temporal  rela¬ 
tions  of  input  patterns  into  spatiotemporal  dynamics  of 
the  output  activity.  We  follow  the  idea  proposed  by 
Hopfield  (1995)  and  assume  that  the  precise  time  ad¬ 
vance  of  a  pattern’s  firing  encodes  the  concentration  of 
an  odor,  though  it  is  a  simplification  of  the  real  olfactory 
code  that  is  yet  to  be  discovered.  However,  the  idea  of 
the  temporal -to-spatial  mapping  could  still  be  applicable 
if  the  temporal  structure  carried  some  other  functional 
significance. 

We  assume  that  spatiotemporal  patterns  in  the  ol¬ 
factory  bulb  are  formed  in  the  following  way:  the  greater 
the  concentration  of  the  odor  component  applied,  the 
earlier  the  correspondent  neural  ensemble  synchronizes 
its  activity  and  fires  (Hoshino  et  al.  1998),  and  the 
greater  is  the  time  q>j  that  the  jth  ensemble  fires  in  ad¬ 
vance  of  the  moment  of  the  maximum  of  its  subthresh- 
old  activation,  which  serves  as  a  reference  time  (Hopfield 
1995), 

According  to  the  Hopfiekfs  (1995)  hypothesis,  the 
corresponding  concentrations  cj  of  n  constituent  mole¬ 
cules  are  encoded  as  time  advances  ^ . . .  t  of  the 
ensemble’s  firing: 

<pj  =  tj- 1{,)  (i) 

where  tj  is  the  time  of  the  ensemble’s  spike  and  /('l  is  the 
reference  time  mentioned  above.  The  functional  relation 
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between  stimulus  intensity  and  time  advance  of  the 
spikes  has  been  proposed  by  Hopfield: 

<Pj  =  *ln(c;/<5)  (2) 

where  each  lime  advance  is  assumed  to  be  proportional 
to  the  logarithm  of  the  corresponding  concentration,  cc  is 
a  coefficient,  and  3  is  a  scale  factor  (Hopfield  1 995). 

Such  logarithmic  scaling  makes  the  relative  time  ad¬ 
vances  of  spatial  patterns  invariant  to  different  concen¬ 
trations  of  the  same  odor.  The  changing  of  the 
concentration  of  a  multicomponent  odor  results  in  a 
lime  advance  of  the  whole  pattern,  while  the  relative 
time  advances  remain  constant. 

2.2  Olfactory  cortex 

2.2.1  Synaptic  organization.  The  piriform  cortex  (PC), 
the  part  of  the  olfactory  cortex  we  focus  on,  is  divided  in 
three  functional  areas;  ventral  and  dorsal  parts  of  the 
anterior  PC  (APCv  and  APCd),  and  the  posterior  PC 
(PPC). 

The  input  (spatial  activation  patterns)  from  the  OB  is 
delivered  directly  to  the  APCv  by  the  lateral  olfactory 
tract  (LOT),  and  via  LOT  collaterals  further  to  the  PPC. 
The  conduction  velocities  along  the  LOT  (7.0  m/s)  are 
greater  than  along  its  collaterals  (1.6  m/s)  (Wilson  and 
Bower  1992).  Thus,  there  is  a  time  delay  of  about  5-7  ms 
between  the  arrivals  of  the  OB  signal  to  anterior  and 
posterior  parts  of  the  PC,  as  measured  by  Ketchum  and 
Haberly  (I993a,b). 

In  addition  to  afferent  input  from  the  OB,  pyramidal 
cells  also  receive  association  projections  from  each  oth¬ 
er,  which  are  distinguished  by  the  area  they  come  from: 
APCv,  APCd,  and  PPC.  The  striking  feature  of  the  PC  is 
the  spatiotemporal  organization  of  its  afferent  and  as¬ 
sociation  connections:  four  major  fiber  systems,  afferent 
input  from  OB,  and  association  projections  from  APCv, 
APCd,  and  PPC  make  synapses  in  distinct  sublayers  of 
the  PC. 

The  closer  the  source  of  the  signal  is  to  the  LOT, 
the  further  from  the  cell  body  corresponding  fibers 
synapse  and  the  greater  is  the  time  of  the  signal's 
propagation  along  the  dendritic  tree  to  the  cell  body. 
Afferent  axons  synapse  on  the  dendrites  of  pyramidal 
ceils  in  layer  la  (the  closest  one  to  the  cortical  surface). 
Association  axons  from  the  APCv  and  the  APCd 
synapse,  respectively,  in  superficial  (sup  lb)  and  deep 
(deep  lb)  parts  of  layer  lb.  Association  axons  from  the 
PPC  excite  the  same  dendrites  deeper  in  layer  deep  lb 
and  in  layer  III. 

2.2.2  Spatiotemporal  dynamics.  The  temporal  dynamics 
of  the  dendritic  trees  in  response  to  LOT  activation  is 
also  strictly  structured.  It  includes  four  peaks  of  inward 
currents:  excitatory  postsynaptic  current  (EPSC)  in  layer 
la  (caused  by  the  signals  of  afferent  fibers),  disynaptic 
EPSC  in  layer  sup  lb  (mediated  by  association  fibers 
from  the  APCv),  small  EPSC  in  deep  lb  (which  is  due,  in 
part,  to  the  APCd’s  signals),  and  inhibitory  PSC  in  layer 


The  structure  of  this  temporal  sequence  of  the  den¬ 
dritic  inward  currents  proved  to  be  crucial  for  their  in¬ 
tegration  at  the  cell  body.  In  order  to  produce  maximum 
activity,  the  cel!  body  needs  not  only  for  both  afferent 
and  association  signals  to  be  applied,  bul  also  for  them 
to  be  temporally  correlated  in  a  specific  way  (Ketchum 
and  Haberly  J993a,c).  This  idea  is  also  supported  by  the 
results  or  Wilson  et  al,  (1992)  (see  Sect.  6.1).  We  explore 
this  integration  mechanism  in  our  model  and  suggest 
that  it  could  be  employed  by  the  olfactory  cortex  for  the 
recognition  of  complex  odors. 

Another  question  we  focus  on  in  this  paper  is:  how  is 
the  concentration  of  odor  components  encoded  in  the 
cortex?  The  experimental  data  of  Duchamp- Viret  et  al. 
(1996)  show  that  in  a  frog  olfactory  cortex  there  is  a 
special  class  of  neurons  which  do  not  discriminate  well 
between  different  odors,  but  instead  seem  to  encode 
odor  concentration.  The  latency  of  their  bursting  was 
found  to  be  correlated  to  the  odor  concentration:  the 
greater  the  concentration,  the  smaller  is  the  latency.  This 
data  is  in  accordance  with  the  idea  of  Hopfield  (1995), 
and  suggests  that  odor  components  and  odor  concen¬ 
trations  may  be  encoded  by  different  populations  of 
neurons. 

However,  such  “concentration  neurons"  were  not 
found  in  other  olfactory  systems,  where  intensity  en¬ 
coding  may  be  different.  In  various  olfactory  systems 
odor  concentrations  influence  the  temporal  structure  of 
OB  activity  patterns.  When  they  are  injected  to  the  PC, 
the  resulting  temporal  patterns  of  incoming  signals  to 
different  PC  regions  invoke  corresponding  sequences  of 
inward  dendritic  currents  at  the  dendritic  trees,  As  well 
as  the  integration  of  EPSCs  in  the  cell  body  depending 
crucially  on  this  pattern  of  inward  currents  (Haberly 
1998;  Ketchum  and  Haberly  1993b,c),  we  conclude  that 
the  dynamics  of  association-fiber  signals  could  be  cor¬ 
related  with  the  concentrations  of  odor  components.  We 
explore  these  ideas  in  our  model,  where  there  is  a  special 
type  of  neurons  sensitive  to  different  concentrations. 


3  The  model 

Let  an  odor  be  a  mixture  of  two  components  A  and  B. 
For  recognition  of  this  mixture  as  a  whole,  some  kind  of 
AND  logic  gate  ODOR=  A  AND  B  has  to  be  realized 
at  some  level  of  a  recognition  system:  in  our  model  this 
AND  function  is  performed  by  integration  of  afferent 
and  association  signals  by  an  integrate-  and-fire  neuron, 
which  fires  if  these  signals  arrive  within  a  close-enough 
time  window. 

We  define  an  odor  as  a  cluster  of  odor  patterns  that 
have  identical  components  with  different  relative  con¬ 
centrations.  An  odor  recognition  system  has  to  be  rough 
enough  to  be  able  to  recognize  that  different  patterns 
may  belong  to  the  same  cluster.  On  the  other  hand,  it 
has  to  be  sensitive  enough  to  distinguish  slightly  differ¬ 
ent  odors  within  a  duster. 

In  our  model  the  recognition  of  an  odor  is  repre¬ 
sented  by  the  firing  of  the  neurons  of  a  specific  ensemble 
in  a  specific  sequence.  Cluster  recognition  and  fine 
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recognition  are  represented  by  activation  of  different 
neural  subensembles. 

Our  network  consists  of  a  layer  of  leaky  integrate- 
and-fire  neurons  u ,,  interconnected  via  arrays  of  inter¬ 
mediate-delay  neurons  d f  (Fig.  I).  The  neurons  Uj  are 
also  connected  with  the  layer  of  temporal  inputs  jy(f). 
These  inputs  simulate  activity  of  the  OB  and  the  neural 
layer  functionally  corresponds  to  the  olfactory  cortex 
that  receives  and  processes  those  patterns. 

The  periodic  inputs  jy(f),  j  =  represent  n 

spatial  patterns  which  correspond  to  n  components  of 
odor  concentration  vector  C=  {ci,C2,...,c„}.  The  in¬ 
puts  are  presented  as  follows: 

sj(i)  =  +  <Pj  -  kT)  (3) 

where  <5{f)  is  Dirac  delta  function,  T  is  the  signal's 
oscillation  period,  s  is  a  spike's  amplitude,  and  the  time 
advances  <Pj  encode  concentrations  of  constituent  mol¬ 
ecules  according  to  (2).  An  example  of  input  pattern  is 
shown  in  Fig.  1. 

There  are  three  types  of  neurons  in  the  layer.  Each 
neuron  is  characterized  by  its  state,  that  is  the  neuron's 
membrane  potential:  «y(f)  for  the  principal  neurons, 
df(f)  for  the  delay  neurons,  and  Vt'(f)  for  the  selective 
neurons.  The  neurons  and  inputs  are  connected  with 
weights  w(neurl,  w(ddl,  and  w(supi  (Fig.  1). 

As  described  below  in  (4),  a  neuron  uj  receives  cor¬ 
responding  input  signal  Sj{t)  from  the  jth  input  and 
lateral  signals  /£(f)  =  i[rf£'(f))  (operator  £  is  defined 
below)  from  the  activated  neuron  which  are  propa¬ 
gated  and  delayed  by  the  delay  neurons  <%. 

The  operator  £  used  above  maps  the  functions  of  a 
neuron  membrane  potential  to  the  function  of  the  spikes 
/(r)  which  this  neuron  produces.  So,  for  example, 
£[«/(0]  =  lA‘)  where  lj(‘)  *s  equal  to 


and  where  f*hresh  are  the  times  when  the  value  of  mem¬ 
brane  potential  uj  reaches  its  threshold. 

If  a  neuron  receives  a  spike  at  time  („  its  potential 
Uj(t)  is  increased  by  the  weighted  value 

lv(neu0  j Sj[t)dl  —  w^neut^s  , 

if  it  is  a  spike  from  the  input  level,  or 

W<0eur>  /  /i'(f)df  =  w<nn,r>  /  £[<(/)]  df  =  w<"eur>S  , 

if  it  is  a  spike  propagated  though  a  delay  neuron  rf(f.  If 
the  potential  of  a  neuron  reaches  its  threshold  value 
u thresh ■  the  neuron  fires.  Its  output  signal  lj{t)  =  £[«,-(<)) 
produces  a  spike  which  is  propagated  to  the  array  of 
delayneurons  that  transfer  the  spike  to  all  other  neurons 
in  the  layer.  At  the  same  time  its  potential  uj  is  instantly 
reset  to  0,  as  shown  in  (5).  Additionally,  the  potential  Uj 
is  constantly  decreasing  with  decay  coefficient  k.  These 
mechanisms  are  employed  by  all  the  neurons  in  the 
model: 

^  =  -*«,{')  +  "(neuS<o  +  ££^’^[<(01(4) 

bi+J  *=l 

uiU  )  =  Mlhr«sh  ^  ui{*+)  =  0  (5) 

The  parameters  of  the  equation  are  set  in  such  a  way 
that  in  order  for  a  neuron  uj  to  fire  it  needs  to  receive 
two  spikes  in  the  narrow  time  window  one  spike, 
Sj(l),  from  the  corresponding  input;  and  another, 
/(f)  =  £[d*  (0]»  from  one  of  the  delay  neurons  (see 
details  in  Sect.  5).  An  exception  is  made  for  the  very  first 
input  spike  in  the  first  cycle,  which  alone  is  able  to 
activate  the  corresponding  neuron.  This  adds  to  the 
model  the  functional  property  of  the  networks  like 
the  so-called  LEGION,  where  the  global  inhibition  of 
the  neurons  depends  on  the  number  of  the  activated 
neurons  (Campbell  and  Wang  1998).  Such  inhibition 
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Fig.  I.  Network  architecture.  The  neurons  of 
activated  sub-ensemble  are  shown  in 

bold *  Arrows  and  black  circles  represen  I  excit¬ 
atory  and  inhibitory  connections,  respectively 
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ensures  that  when  fewer  neurons  are  activated,  the 
greater  is  the  probability  for  a  neuron  to  fire.  In  our 
model  the  neurons  are  made  more  sensitive  to  the  very 
first  input  spike,  because  there  is  no  activation  yet  in  the 
neuron  layer.  The  biological  correspondence  of  this 
assumption  is  discussed  in  Sect.  6.1. 

Delay  neurons  in  the  arrays  are  integrate-and-fire 
neurons  with  added  inherent  propagation  delays  D 
k  =  1,...  ,n,  defined  as 

=  -*<(/)  +  w<wl£[«,(f  -  D,)| 

+  £  wl“w£|*f</)|  (6) 

The  parameters  of  the  equation  make  the  delay  neuron 
work  as  a  delay  operator  (details  in  Sect.  5).  When 
receives  a  spike  at  time  t  it  fires  at  r  +  Dk.  Functionally, 
dff  are  equivalent  to  the  ordinary  integrate-and-fire 
neurons  as  in  (4),  connected  to  the  principal  neuron  ut 
via  corresponding  delay  The  values  of  delays  Dk 
change  gradually  across  the  array  as  follows: 


Dk  = 


T(k-\/l) 

m  * 


k  =  l  t  * . .  r  m 


(7) 


where  T  is  the  oscillation  period  (3),  and  m  is  the  number 
of  the  delay  neurons  in  the  array.  This  specific  delay 
distribution  ensures  recognition  properties  of  the  system 
which  are  discussed  below. 

As  a  pattern  processing  example  we  consider  a  simple 
case  where  an  odor  with  two  components 
{0, . . . , Ci,Cj,..  . , 0}  is  applied  with  c/  >  cj.  For  conve¬ 
nience  we  will  be  presenting  this  as  {c/,cj}.  When  this 
odor  is  applied,  the  neuron  «/  receives  the  input  spike 
s/(t)  first,  at  the  moment  ft,  and  then  uj  receives  sj(t)  at 
h- 

The  single  spike  s/  is  enough  for  u/  to  fire  because  of 
the  exception  mentioned  above.  The  neuron  fires  at 
the  moment  ft  and  sends  spikes  I[w/(/)]  to  all  other 
neurons  «/>  j  —  1,. . . ,  n,  j  ^  1  via  delay  arrays  d f.  The 
neurons  us  receive  the  delayed  signals  from  df  at  times 
/  +  £>*,  k  =  ,m.  However,  this  is  not  enough  for 

them  to  fire  because  a  spike  from  the  input  level  is  also 
needed,  Although  all  of  them  show  subthreshold  acti¬ 
vation,  only  the  neuron  uj  which  will  receive  the  input 
spike  sj  will  actually  reach  the  threshold  and  fire.  Fi¬ 
nally,  the  neurons  that  fire  are  u/t  uj,  and  all  interme¬ 
diate  neurons  df  in  the  array  which  connects  them. 

The  values  of  £>*  set  by  (7)  ensure  that  one  and  only 
one  of  the  delay  neurons  fires  and  sends  a  spike  to  the 
neuron  uj  within  the  time  window  A/w.  Thus,  although 
all  of  the  delay  neurons  in  the  array  fired,  only  one  of 
them  actually  contributes  to  the  firing  of  the  neuron  uj. 
To  distinguish  this  contributing  delay  neuron  from 
others,  an  additional  layer  of  selective  neurons  jc£  is 
added  (Fig.  1). 

Selective  neurons  jr*  are  functionally  equivalent  to  the 
principal  neurons  uj.  They  work  as  coincident  time  de¬ 
tectors  and  fire  if  spikes  from  neurons  uj  and  df  arrive  at 

within  time  window  Ac 


=  -fcq(r)  + 

+  w(neur>I[«,(f)3  (8) 

Selective  neurons  provide  negative  feedback 
to  the  delay  neurons  df,]  that  did 
not  contribute  to  the  firing  of  uj  (6).  Because  of  this 
selective  feedback,  neurons  df  which  contributed  to  the 
firing  of  uj  stay  unchanged,  while  the  rest  of  the  delay 
neurons  are  suppressed  and  will  not  be  sensitive  to  the 
spikes  from  during  suppression  time  7s,  defined  as 
follows: 


Ts 


1  /  w(*up>j  \ 


(9) 


At  the  output  level  there  is  now  a  sequence  of  firing  of 
neurons  «/,  df[,  and  «/  one  after  another.  Firing  of  the 
neurons  «/  and  uj  indicates  that  an  odor  of  the  cluster 
{c/,o}  is  recognized,  The  firing  of  the  specific  delay 
neuron  df  defines  the  relative  concentration  of  two 
components,  the  logarithm  of  which  lays  in  the  vicinity 
of  the  delay  DK  of  the  corresponding  delay  neuron: 


DK  -  —  <  <x\n 


£/ 

CjJ 


<  Dk  +  t— 
2m 


(10) 


An  example  of  the  system’s  dynamics  is  shown  in  Figs.  2 
and  3.  During  the  first  cycle  (0  <  t  <  T),  neurons  df  are 
activated  by  «/  and  fire  consequently  in  the  order  of  their 
inherent  delays  Dk,  starting  from  the  one  with  the 
smallest  Dk,  until  df  gets  activated,  which  fires  within 
time  window  At  with  uj.  The  correlated  firing  of  uj  and 
df  makes  xf  fire.  Selective  neuron  xf  suppresses  the 
remaining  delay  neurons  for  the  suppression  time  7s, 
which  is  equal  to  the  oscillation  period  T  =  20.  After 
time  T,  selective  neuron  xf  will  suppress  them  again.  So, 
the  only  delay  neuron  that  will  fire  during  the  second 
and  following  cycles  is  df.  The  parameters  used  in  the 
simulation  are  described  in  Sect.  5. 


Fig.  2-  Dynamics  of  the  neuraJ  ensemble  during  the  recognition 
process.  The  bars  represent  the  spikes  of  the  corresponding  neurons 
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Fig.  3a, b.  Simulation  of  spatioiemporal  dynam¬ 
ics  of  ihc  ensemble;  a  distribution  of  the  neurons 
in  the  layer:  b  spa tio temporal  neural  dynamics 
Axis  z  corresponds  to  the  membrane  potentials 

w(0  or  d(r)-  k  =  1 , 2 . N;  and  7*  is  the  period  of 

the  oscillations 


4  Temporal  segmentation 

According  to  Hopfield’s  (1985)  mechanism  of  formation 
of  the  input  temporal  sequence,  if  a  mixture  of  several 
odors  is  applied,  the  corresponding  spatiotemporal 
sequences  are  superimposed  and  the  resulting  input 
sequence  contains  patterns  of  all  components  of  each 
odor.  The  stronger  the  component,  the  earlier  its  pattern 
fires,  no  matter  which  odor  the  component  belongs  to* 
In  order  to  segregate  the  odor  patterns  in  time,  one  of 
the  ensembles  should  win  and  hence  suppress  the  others 
for  a  period  of  several  cycles.  After  that,  due  to  the 
neural  fatigue,  the  winning  ensemble  stops  firing  and  the 
second- strongest  ensemble  wins  and  fires  during  the  next 
several  cycles  (Campbell  and  Wang  1998;  Hoshino  et  ah 
1998;  Malsburg  1992).  To  realize  such  pattern  segmen¬ 
tation,  the  modified  network  from  Sect.  2  with  addi¬ 
tional  neural  interaction  and  neural  fatigue  function 
F{p)  is  used: 

+EE^"r)iKW] 

+  £  wjin,er)i[Wi(/)]}F(p)  tio 

rXi.UiUUj&E, 

where  F(p)  is  defined  as  a  step  function,  with  F(p)  —  1  if 
2kpf  <p<(2k+\)pf,  and  F  (p)  =  0  if  (2k+\)pf 
<  p  <  (2k  +  2  )pjr\  VA  :  k  =  1,2  The  variable  p  is 

the  number  of  times  the  potential  of  neuron  uj  reached 
its  threshold,  and  pr  is  the  number  of  firings  after  which 
the  neuron  becomes  insensitive  to  the  inputs  and  stays 
silent  for  another  p/  cycles. 


Neural  ensemble  Et,  z  =  I, . . .  ,Z,  is  a  group  of  neu¬ 
rons  that  correspond  to  the  components  of  the  same 
odor.  All  neurons  which  are  not  members  of  the  same 
ensemble  are  interconnected  with  negative  weights 
timer) ,  ^yhen  a  neuron  «,■  fires,  it  sends  inhibitory  signals 

[tii(01  to  [he  neurons  that  do  not  belong  to  any  of 
the  ensembles  E.  in  which  the  neuron  ut  participates.  The 
membrane  potentials  of  those  neurons  are  decreased  by 
wj;n'“>£,[u,(0].  They  stay  suppressed  for  the  refractory 
period  7r  during  which  they  cannot  fire  regardless  of  the 
input  signals  received.  7r  is  determined  as  follows: 


7k 


=M 


wtipter>J  \ 

MlhtBh  - 


(12) 


The  dynamics  of  the  competition  between  ensembles  is 
quite  a  complicated  process  (Campbell  and  Wang  1998; 
Hoshino  et  ai.  1998),  because  Tor  a  neuron  in  the 
ensemble,  its  probability  of  being  suppressed  depends 
on  the  statistical  value  of  the  difference  of  received 
inhibitory  and  excitatory  spikes. 

In  our  model,  according  to  (1 1)  and  the  exception  for 
the  first  spike  in  the  input  sequence  (Sect.  3),  the  first 
activated  neuron  suppresses  the  others  which  are  not  in 
its  ensemble,  and  does  not  allow  them  to  fire  with  100% 
probability.  So,  the  ensemble  which  contains  the  win¬ 
ning  neuron,  or,  in  other  words,  the  odor  with  the 
strongest  component,  always  wins  the  competition  first. 

As  an  example  we  consider  the  case  where  two  odors 
{c/,o}  and  {c/>,cg}  are  applied  to  the  network  with  the 
following  order  of  the  concentrations  of  their  compo¬ 
nents:  ci  >  cp  >  cj  >  eg.  The  temporal  sequence  of  input 
spikes  is  the  same,  {st(t),sp(t),Sj(t),SQ(t)},  as  shown  in 
Fig.  4.  The  first  of  them,  s/(t),  activates  neuron  w/,  which 
sends  inhibitory  signals  to  up  and  uq,  and  excitatory  signal 
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to  uj  via  an  array  of  the  delay  neurons.  When  input  spike 
jp{f)  appears,  the  corresponding  output  neuron  up  is 
suppressed  and  will  not  respond  to  the  impulse.  Then  in¬ 
put  sj  activates  neuron  uj,  which  already  received  an  ex¬ 
citatory  signal  from  u/;  After  that  sq(i)  fails  to  activate 
neuron  uq.  Finally  the  neurons  u ,  and  uj  fire,  while  uP  and 
uq  remain  silent,  In  this  way  the  odor  {w/,u/}  is  segmented 
from  the  background  and  attention  is  focused  on  it  for  the 
period  of  pr  =  2  cycles  (Fig,  3).  Then,  due  to  the  neural 
fatigue  F\p\,  the  neurons  w/  and  uj  stop  firing  and  up  and 
uq  become  activated,  so  the  attention  is  now  refocused  on 
this  odor.  So,  the  odor  patterns  are  temporally  segregated 
and  processed  one  at  a  time,  as  is  shown  in  Figs.  4  and  5. 
The  parameters  used  in  the  simulation  are  specified  in 
Sect.  5. 


5  Simulation 

The  neurons  described  by  (4)  and  (8)  work  as  coincident 
time  detectors.  They  fire  if  two  spikes  arrive  at  a  neuron 
within  time  window  At,  the  size  of  which  is  defined  by 
the  parameters  of  the  integrate-and-fire  neurons  as 
follows: 

k  \w^)s  )  j 

k  vy11*"1-*  ) 

where  A/to  and  A/to  are  the  time  windows  of  neurons  «; 
and  jcf ,  respectively.  In  our  simulation  uthr«sh  =  xihresh.  so 
A(i“)  — -  A/to,  and  we  will  represent  them  both  as  At. 

Since  this  is  essentially  the  property  of  the  neurons 
used  in  the  model,  there  was  no  need  to  implement  them 
by  the  actual  integrate-and-fire  neurons.  The  neurons  u; 
were  replaced  by  logic  units  uj,  characterized  by  its  state 
u}(/)  that  has  the  following  properties.  Without  inputs, 


k*(0  is  equal  to  0.  If  uj  receives  two  spikes  with  weights 
wtne“0  at  the  moments  f]  and  with  ti  >  /t,  then 

if  <  At 

then  «’(£)  =  «,hmh 

else  u*  (/)  =  0  (14) 

If  the  state  uj(r)  of  logic  unit  uj  reaches  its  threshhold 
value  ujhresh  at  time  /thresh,  it  behaves  in  the  same  way  as 
neurons  u.  do.  The  value  of  u*{/)  is  reset  to  0  and  the 
output  function  of  the  unit,  £[«}(/)],  is  equal  to 

A(/  -  /,hreih)-  The  rules  for  units  x'kJ‘  are  analogous  to 
the  rules  of  uj.  The  state  of  cT/  is  defined  as 

The  focus  of  the  mode)  is  the  computational  abilities 
of  spatiotemporal  integration  and  temporal  dynamics 
discussed  in  Sect.  3.  For  this  reason,  temporal  parameters 
have  been  assigned  biologically  realistic  values,  such  as 
the  period  of  oscillations,  T  —  20  ms,  which  corresponds 
to  50- Hz  oscillations  observed  in  the  olfactory  cortex, 
This  also  makes  the  values  of  time  delays  DK,  defined  by 
(7),  comparable  with  the  delays  of  5-7  ms,  related  to 
different  association  fibers  (see  Sect.  6. 1  for  details). 

The  values  of  the  following  parameters  were  chosen 
arbitrarily:  T  =  Tr  =  7s  =  20,  A/  =  5,  a  =  4,  5=1, 

s  =  1 »  uihresh  =  *ihresh  —  thresh  =  1 ,  Pf  —  2,  n  =  4,  and 
m  =  4. 

The  weights  w(del)  =  1.1  were  assigned  their  values  in 
order  to  make  a  single  weighted  spike  nv(dt,)  be  enough 
to  fire  a  delay  neuron.  The  values  of  weights 
^(neur)  _  q  75  wcre  defined  so  as  to  ensure  that  not  one, 
but  two  weighted  spikes  sw(neur)  (and  in  a  small-enough 
time  window)  are  needed  to  activate  a  principal  neuron. 
The  parameters  w<sup)  =  -8.0963,  =  -39.91,  and 

k  —  0.2197  were  defined  by  (9),  (12),  and  (13),  respec¬ 
tively. 

The  dynamics  of  an  example  simulation  is  shown  in 
Fig.  5,  where  the  mixture  of  odors  {u/,i//}  =  {100,50} 


Fig.  4.  Dynamics  or  two  neural  ensembles  during 
processes  of  temporal  binding,  segmentation,  and 
attention 


65 


a 


t  =  S  ^  kT 


1  =  1 9  +<  k+prJT 


Fig.  5a,b.  Simulation  of  the  neural  dynamics  of 
pattern  segmentation  and  attention  focusing: 
a  distribution  of  the  neurons  in  the  layer; 
b  spa ti ©temporal  neural  dynamics  during  the  first 
p/  ( left  column)  and  second  p/  {right  column) 
cycles.  Axis  z  corresponds  to  the  membrane 
potentials  w(r)  or  k  =  1(. . ,  ,pf  t  and  T  is  the 
period  of  the  oscillation 


and  {up,uq}  =  {S0T  3}  is  applied.  During  the  first  three 
cycles  the  neurons  of  the  first  odor  fire  in  the  following 
sequence:  and  uj  (Fig.  5b,  left  column).  This  in¬ 

dicates  that  the  odor  {w/Tt//}  is  recognized  with  the 
relative  concentrations  of  its  components  defined  by  (9) 
as  0  <  In—]  <  1.25,  After  three  cycles  {up.d^1  ,uq}  be¬ 
comes  activated  and  suppresses  in  its  turn  the  first  en¬ 
semble  (Fig.  5b,  right  column).  So,  the  ratio  of  the 
concentrations  of  the  components  of  this  odor  is 
2.5  <  In^  <  3.75. 

6  Discussion 

6,1  Biological  correspondence 

We  argue  that  the  dynamics  of  our  model  reflects  certain 
aspects  of  olfactory  cortex  activity  related  to  informa¬ 
tion  processing.  The  key  principles  of  mixture  recogni¬ 
tion  in  our  system  are  that  a  single  input  spike  from  the 
OB  is  enough  to  activate  a  certain  neuron  (e.g.,  uh  in  the 
model)  while  the  others  (e.g.,  uj)  need  signals  both  from 
the  OB  and  from  another  principal  neuron,  propagated 
by  the  delay  neurons. 

Experimental  results  (Wilson  and  Bower  1992)  indeed 
suggest  that  while  anterior  areas  of  the  cortex  may  be 
activated  by  afferent  input  only,  the  posterior  parts  re¬ 
quire  both  afferent  and  association  signals.  This  mech¬ 
anism  rises  from  the  properties  of  the  relative  density  of 
association  and  afferent  terminals,  defined  by  the  asso- 
riation-to-afferent  dominance  ratio  (AADR).  The 
A  ADR  increases  along  the  distance  from  the  LOT:  in 
the  APCv  (the  closest  area  to  the  LOT),  afferent  fibers 


prevail  over  the  association  ones,  and  the  opposite  is 
true  for  the  PPC  (Haberly  1998).  The  results  of  Wilson 
and  Bower  (1992)  show  that  when  the  afferent  activity  is 
not  strong  enough,  it  can  only  activate  the  APC,  where 
the  AADR  is  low.  Activation  of  the  posterior  cortex 
(with  a  high  AADR)  results  mostly  from  delayed  arrival 
of  association- fiber  activity.  These  results  suggest  that 
the  areas  with  medium  AADRs,  such  as  the  APCd.  may 
need  equally  both  afferent  and  association  inputs  to  be 
activated. 

In  Sect.  2.2.2  we  described  the  spariotemporal  se¬ 
quence  of  the  dendritic  currents,  produced  by  afferent 
and  association  signals.  Its  correspondence  with  the 
dynamics  of  our  model  is  the  following:  neuron  u;  re¬ 
ceives  the  spikes  Sj(t)  from  the  OB,  and  L[d£]  from  the 
principal  neurons,  differently  delayed  by  delay  neurons 
(or  association  fibers,  as  in  the  real  PC),  as  shown  in 
Fig.  2.  Each  of  the  incoming  signals  increases  the  state 
of  the  neuron  (i.e.,  produces  a  peak  of  inward  current). 
Uj  is  activated  only  if  afferent  input  sj{t),  and  one  of  the 
association  inputs  L[d£ j,  arrive  in  a  small-enough  time 
window  (in  terms  of  the  PC,  the  induced  EPSCs  are 
close  enough  in  time). 

This  mechanism  is  supported  by  experimental  results 
described  in  Ketchum  and  Haberly  (!993b,c)  and  Ha¬ 
berly  (1998),  which  show  that  the  optimal  temporal 
correlation  of  inward  currents  is  the  one  which  ensures 
their  roughly  synchronous  arrival  at  the  cell  body.  In 
these  experiments,  when  three  EPSCs  were  presented  in 
their  natural  sequence  (with  interpeak  latencies  of  sev¬ 
eral  milliseconds),  the  induced  potential  at  the  cell  body 
was  50%  greater  than  the  one  caused  by  the  same 
EPSCs  applied  simultaneously. 
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Moreover,  the  delays  produced  by  propagation  along 
the  dendritic  tree  could  synchronize  the  incoming  sig¬ 
nals,  as  did  delay  neurons  in  our  model.  Indeed,  different 
fibers  synapse  at  distinct  distances  from  the  ceil  body, 
and  thus  require  different  times  for  a  postsynaptic  spike 
to  propagate  to  it  via  the  dendritic  tree.  The  latency 
between  peak  depolarization  in  afferent  input  in  layer  la 
and  the  peak  depolarization  of  the  cell  body  is  approx¬ 
imately  6  ms  (Haberly  1998). 

The  repetitive  temporal  dynamics  in  our  model  is  also 
related  to  the  one  of  the  real  PC,  where  the  temporal 
sequence  of  dendritic  postsynaptic  currents  evolves  in 
each  cycle  of  50-Hz  cortical  oscillations.  Odor  compo¬ 
nents  are  also  segregated  temporally,  with  each  of  them 
processed  during  different  oscillation  cycles.  Thus,  the 
temporal  sequence  of  EPSCs,  whatever  functional  role  it 
plays,  participates  repetitively  in  the  processing  of  each 
odor  component. 

Selective  suppression  of  the  delay  neurons  in  our 
model,  which  corresponds  to  the  suppression  of  associ¬ 
ation-fiber  signals,  also  has  its  prototype  in  the  PC.  As 
shown  in  Ketchum  and  Haberly  (1993a),  when  two 
successive  shocks  are  applied  to  the  LOT,  the  first  one 
induces  the  full  sequence  of  the  peaks  of  inward  cur¬ 
rents,  but  the  second  one  produces  only  an  isolated 
monosynaptic  EPSC  coming  from  the  APCv.  The  as¬ 
sociation-fiber  signals  are  blocked  by  inhibitory  postsy- 
naptic  currents  evoked  by  the  first  shock.  However, 
there  is  not  enough  biological  evidence  yet  to  specify 
functional  significance  of  this  suppression  and  compare 
it  to  the  one  of  our  model,  where  all  delay  neurons  are 
active  at  the  first  cycle  and  get  selectively  suppressed  at 
the  second  one  (Fig.  2). 

In  our  model,  activity  of  a  delay  neuron  represents 
corresponding  concentration  of  the  odor  component. 
Such  behavior,  although  detected  in  the  frog  olfactory 
cortex  (Duchamp- Vi  ret  et  ai.  1996),  was  not  observed  in 
other  olfactory  systems.  We  argue  that  the  dynamics  of 
delay  neuron  in  our  model  can  be  seen  as  differently 
delayed  activity  of  association  fibers  coming  from  dif¬ 
ferent  cortical  areas:  APCv,  APCd,  and  PPC.  They  also 
are  of  the  same  range:  the  time  of  signal  propagation 
along  the  cortex  -  from  anterior  to  posterior  parts  -  is 
5-7  ms. 


6.2  Olfactory  system  models 

The  general  idea  behind  most  of  the  odor  recognition 
models  is  that  an  odor  pattern  is  temporally  segregated 
in  the  patterns  of  its  components,  which  are  processed  in 
the  order  of  their  significance  or  intensity  (the  strongest 
one  is  processed  first).  However,  it  is  not  absolutely  clear 
which  parts  of  the  olfactory  system  are  involved  in  this 
process,  and  to  what  degree. 

In  the  model  of  Ambros-Ingerson  et  al.  (1990),  tem¬ 
poral  segmentation  of  the  bulbar  activation  patterns  is 
realized  by  the  selective  inhibitory  cortical  feedback. 
During  each  cycle,  OB  activity  is  projected  to  the  cortex 
where  the  pattern  of  its  strongest  component  suppress 
the  other  ones  (the  winner-takes-all  mechanism)  and 


activates  inhibitory  feedback  to  the  bulb.  This  feedback 
suppresses  the  part  of  the  bulbar  activity  which  corre¬ 
sponds  to  this  strongest  component.  During  the  fol¬ 
lowing  cycles,  the  normalized  remainder  of  the  input  is 
presented  to  the  cortex  and  the  same  process  occurs 
repetitively.  During  different  cycles,  input  is  classified  as 
a  part  of  a  cluster  of  the  corresponding  hierarchical 
level.  Processing  during  the  first  cycle  corresponds  to  the 
first  level  clusters  (rough  recognition),  subcluster  recog¬ 
nition  is  realized  at  the  second  cycle,  and  so  on. 

Input  segmentation  in  this  model  is  due  to  the  corti- 
cobulbar  interaction.  However,  it  may  be  also  realized  in 
the  OB  without  corticobulbar  interplay.  In  the  model  of 
Hoshino  et  al.  (1998)  the  winner-takes-all  competition 
occurs  in  the  OB  itself.  The  winning  pattern  suppresses 
the  others,  stays  active  during  several  cycles,  and  then 
stops  firing  due  to  the  neural  fatigue.  The  second- 
strongest  pattern  then  wins,  and  the  process  repeats. 
Chaotic  dynamics,  which  is  believed  to  be  crucial  for 
bulbar  information  processing,  is  also  taken  into  ac¬ 
count  in  this  system. 

The  models  discussed  above  deal  with  odors  which 
have  distinct  components.  If,  instead,  complex  odors 
with  the  same  components  but  in  different  concentra¬ 
tions  are  ideally  mixed,  the  information  about  them  is 
lost,  and  their  separation  is  impossible.  An  odor  mixture 
{4A,  7B}  may  be  composed  of  an  infinite  number  of 
odor  combinations,  such  as  {2A,  B}  +  {2A,  6B}  or  {A, 
3B}+  {3A,  4B}.  However,  the  problem  becomes  solv¬ 
able  if  temporal  fluctuations  of  the  odors  are  indepen¬ 
dent.  This  is  the  case  when,  for  example,  the  sources  of 
odors  are  spatially  distinct  and  there  is  enough  turbu¬ 
lence  in  the  airflow.  This  idea  was  explored  in  the  models 
of  Hendin  et  al.  (1994,  1998),  where  temporal  segmen¬ 
tations  of  the  bulbar  activity  is  based  on  the  temporal 
fluctuations  of  input  odors. 

In  our  model  we  focus  on  the  recognition  of  the  odor 
mixtures  which  are  already  segmented  into  spatiotem- 
poral  patterns  of  OB  activity,  after  the  temporal  fluc¬ 
tuations  of  odors,  if  any,  were  employed  by  the  OB.  We 
propose  a  new  possible  functional  role  of  temporal 
correlation  of  association-fiber  signals,  activity  of  which 
is  influenced  by  the  temporal  structure  of  the  OB  signals, 
which,  in  turn,  is  assumed  to  be  correlated  with  the  odor 
concentrations. 

The  advantage  of  the  distributed  representation  of 
our  model  is  its  flexible  recognition  ability.  Speaking  of 
the  olfactory  system,  where  the  number  of  odors  and 
their  mixtures  the  brain  can  possibly  perceive  is  enor¬ 
mous,  one  of  the  computational  problems  is  how  to 
preserve  the  hierarchy  and  similarity  in  the  represen¬ 
tation  of  odor  memory.  In  our  system  the  neural  en¬ 
semble  which  corresponds  to  the  odor  of  coffee  would 
be  a  part  of  a  bigger  ensemble  that  represents  more 
complicated  odors  of  a  coffee  house.  The  dynamics  of 
the  smaller  ensemble  also  would  be  part  of  the  bigger 
ensemble’s  dynamics.  On  the  other  hand,  two  similar 
odors  (odors  with  some  common  components)  would 
activate  two  similar  neural  ensembles  with  close  spatial 
dynamics.  This  would  happen  because  common  com¬ 
ponents  and  their  relative  concentrations  would 
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activate  the  same  spatial  dynamics  of  the  same  en¬ 
semble,  The  neural  processing  of  spatio temporal  pat¬ 
terns  is  one  of  the  key  principles  for  the  brain's  fast 
flexible  pattern  recognition  that  is  still  beyond  the  reach 
of  modem  computational  techniques.  In  this  paper  we 
present  an  attempt  to  understand  the  principles  of  how 
the  brain  could  use  its  spatiotemporal  dynamics  for 
recognition  tasks. 
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Abstract.  This  paper  presents  a  model  of  a  network  of  integrate- a  nd-fire  neurons  with  time 
delay  weights,  capable  of  invariant  spatio-temporal  pattern  recognition.  Spatio-temporal  pat¬ 
terns  are  formed  by  spikes  according  to  the  encoding  principle  that  the  phase  shifts  of  the 
spikes  encode  the  input  stimulus  intensity  which  corresponds  to  the  concentration  of  con¬ 
stituent  molecules  of  an  odor  We  applied  the  Hopfield’s  phase  shift  encoding  principle  at  the 
output  level  for  spatio-temporal  pattern  recognition:  Firing  of  an  output  neuron  indicates  that 
corresponding  odor  is  recognized  and  phase  shift  of  its  firing  encodes  the  concentration  of  the 
recognized  odor.  The  temporal  structure  of  the  model  provides  the  base  for  the  modeling  of 
higher  level  tasks,  where  temporal  correlation  is  involved,  such  as  feature  binding  and  seg¬ 
mentation,  object  recognition,  etc. 

Key  words,  integrate-and-fire  neurons,  olfactory  cortex.,  phase  shift  encoding,  spatio-temporal 
pattern  recognition 


I-  Introduction 

1. 1.  OLFACTION 

In  the  olfactory  bulb  odor  stimulus  information  is  encoded  into  a  periodic  spatio- 
temporal  pattern  of  oscillatory  neural  activity  {Hoshino  et  af,  1998  Skarda  and 
Freeman,  1987;  Laurent  and  Davidowitz,  1994;  Laurent,  1996),  Specific  spatial  pat¬ 
terns  of  synchronized  firing  are  correlated  to  a  certain  constituent  molecules  of  the 
applied  odor.  Relative  time  advances  of  the  appearances  of  these  spatial  patterns  are 
correlated  with  the  concentrations  of  the  molecules.  This  spatio-temporal  correlation 
enables  odor  recognition  and  concentration  estimation  in  the  olfactory  cortex  (Hop- 
field,  1995;  Hoshino  et  al,,  1998).  The  olfactory  epithelium  of  a  nasal  cavity  contains 
hundreds  (say  n )  of  receptors  sensitive  to  various  types  of  molecules  (Rcssler,  1994). 
Thus*  epithelium  perceives  odors,  basically,  as  the  mixtures  of  their  components 
(different  molecules).  It  is  convenient  for  the  odors  to  be  represented  by  the  con¬ 
centration  vector  C  =  {c{ ,  C2,  * .  - ,  cn}>  where  Cj  is  the  corresponding  concentration  of 
the  jth  odor  component. 

During  sniffing  constituent  molecules  of  an  odor  cause  receptor  neurons  (that  are 
equally  distributed  in  the  olfactory  epithelium)  to  fire,  producing  spikes  (Duchamp- 
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Viret  et  al.,  1998;  Ressler,  1994).  A  firing  pattern  in  the  receptor  field  is  then  mapped 
to  the  olfactory  bulb.  Axons  from  the  same  receptors  of  olfactory  epithelium  go  to 
the  corresponding  ensembles  of  the  neurons  in  the  olfactory  bulb  (Ressler,  1994).  So, 
each  ensemble  of  neurons  in  the  olfactory  bulb  represents  a  corresponding  consti¬ 
tuent  molecule  of  the  odor.  Once  they  receive  impulses  from  the  receptor  level,  the 
neurons  in  the  ensembles  synchronize  their  activity.  They  fire  synchronously  within 
each  ensemble,  and  with  the  time  shifts  between  different  ensembles  (Hoshino  et  al., 
1998;  Laurant  and  Davidowitz,  1994).  As  a  result,  the  odors  at  the  bulbar  level  are 
presented  by  the  sequence  of  firing  of  different  neural  ensembles,  where  the  specific 
neural  ensembles  represent  the  odor  constituent  molecules  and  the  time  shifts  of  the 
firing  of  these  ensembles  represent  the  concentrations  of  the  constituent  molecules 
(Hoshino  et  al.,  1998;  Duchamp- Viret  et  al.,  1996). 

The  larger  is  the  concentration  c/  of  a  constituent  molecule,  the  earlier  the  cor¬ 
responding  ensemble  gets  excited  and  synchronized  (Campbell  andWang,  1998)  and 
the  greater  is  the  time  <pj  the  yth  ensemble  fires  in  advance  of  the  moment  of  the 
maximum  of  its  subthreshold  activation,  which  serves  as  a  reference  time  (Hopfield, 
1995).  So,  the  corresponding  concentrations  c/  of  n  constituent  molecules  can  be 
encoded  as  time  advances  <p,,g>2 . <p„  of  ensemble’s  firing  presented  as: 

<Pj  =  {J~  (1) 

where  tj  is  the  time  of  the  ensemble’s  spike  and  is  the  reference  time  mentioned 
above.  The  relationship  (2)  between  input  concentrations  cj  and  corresponding  time 
advances  <pj  has  been  proposed  in  the  Hopfield’ s  model  (Hopfield,  1995),  where  each 
time  advance  is  proportional  to  the  logarithm  of  the  corresponding  concentration,  a 
is  a  coefficient  and  <5  is  a  scale  factor. 

<#>,  =  <*  in (cj/S)  (2) 

The  logarithmic  scaling  makes  the  relative  phase  pattern  invariant  to  the  different 
concentrations  of  the  same  odor.  In  this  case,  not  the  relative,  but  the  entire  pattern  is 
shifted  when  the  odor  concentration  is  changed. 

1.2.  OLFACTORY  PATTERN  RECOGNITION  MODELING 

In  most  of  the  functional  olfactory  models  reported  in  the  literature  the  authors  focus 
on  the  spatio-temporal  pattern  formation  at  the  olfactory  bulb  level.  The  basic 
principles  of  pattern  formation  presented  above  are  more  or  less  understood 
and  they  are  supported  by  the  biological  experiments  (Duchamp-Viret  et  al., 
1998;  Hoshino  et  al.,  1998;  Laurent  and  Davidowitz,  1994;  Laurent,  1996).  However, 
the  mechanism  by  which  these  spatio-temporal  patterns  (Figure  1(a))  are  recognized 
and  then  processed  at  the  next,  olfactory  cortex  level,  remains,  basically,  unknown. 

One  of  the  plausible  biological  mechanisms  of  spatio-temporal  pattern  recognition 
is  a  system  with  an  appropriate  set  of  delays  stored  in  synaptic  memory  followed  by 
coincidence  time  detectors  that  receive  these  appropriately  delayed  (and  now  syn¬ 
chronized)  signals.  Recognition  of  a  stored  odor  can  be  indicated  by  firing  of  the 
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corresponding  coincidence  time  detector  neuron  (Hopfield,  1995;  Natschlager  and 
Ruf,  1998;  White  et  al.,  1998).  The  question  that  arises  is  how  the  concentration  of 
the  recognized  multi-component  odor  is  represented  at  the  olfactory  cortex? 

Experimental  evidence  indicates  that  time  shift  encoding  of  stimulus  intensity 
occurs  not  only  at  the  bulb  level,  but  at  olfactory  cortex  level  as  well  (Duch- 
amp-Viret  et  at.,  1996,  1998).  This  supports  the  idea  that  time  shift  encoding  is 
also  used  at  the  'output  level’  of  the  olfactory  system,  i.e.  at  the  olfactory  cortex. 
It  also  indicates  that  the  mechanism  of  this  encoding  could  be  similar  to  the  one  at 
the  olfactory  bulb. 

However,  most  of  the  olfactory  pattern  recognition  models  do  not  make  use  of 
temporal  encoding  and  temporal  processing.  In  such  models  the  patterns  to  be 
recognized  are  ordinary  time  independent  vectors  that  represent  certain  odor  qua¬ 
lities.  Those  vectors  are  then  recognized  by  one  of  the  classical  pattern  recognition 
techniques,  which  do  not  have  much  in  common  with  biological  temporal  processing. 

Temporal  encoding  and  temporal  processing  have  only  recently  been  included  in 
the  olfactory  pattern  recognition  modeling.  Significant  progress  in  this  direction  has 
been  made  by  White  et  al.  (1998).  In  their  mode!  vapor  identity  is  encoded  by  the 
spatial  code  across  output  units,  and  vapor  intensity  is  represented  by  response 
latency.  The  system  is  not  only  biologically  relevant  (at  some  extent),  but  also  proved 
to  be  more  effective  than  classical  neural  networks  models.  Its  percentage  of  correctly 
identified  test  patterns  was  higher  than  the  one  of  the  feed-forward  neural  network 
with  hidden  layer  (82%  and  71%  correspondingly)  (White  et  al.,  1998). 

However,  the  model  of  White  et  al.  is  not  complete.  Odor  intensities  are  encoded 
just  qualitatively:  shorter  response  latency  signifies  greater  concentration  (and  vice 
versa),  but  no  precise  functional  correspondence  exists  between  response  latency  and 
odor  concentration.  Moreover,  the  model  as  it  is  shown  in  (White  et  al.,  1998)  works 
only  in  the  very  narrow  range  of  relative  concentrations.  Our  proposed  model  imple¬ 
ments  precise  functional  encoding  of  the  pattern  intensity,  that  is  odor  concentration, 
with  the  phase  shifts  of  the  output  neuron  firing.  The  odor  recognition  remains 
invariant  within  broad  range  of  concentrations  due  to  the  Hopfield’s  logarithmic 
intensity  encoding. 

2.  Model 

Our  network  consists  of  one  layer  of  m  leaky  integrate-and-fire  neurons  fully  con¬ 
nected  with  n  temporal  inputs.  These  inputs  simulate  spatio-temporal  patterns 
formed  in  the  olfactory  bulb  (Figure  1(a)),  and  the  neural  layer  (Figure  1(b))  cor¬ 
responds  to  the  olfactory  cortex  that  receives  and  recognizes  those  patterns.  The 
periodic  inputs  jr/r)  for  j  =  n  are  expressed  by  Dirac  delta  function: 


(3) 


*=i 


where  T  is  the  period  of  the  signal’s  oscillation.  The  time  advance  <p,-  of  a  periodical 
input  spike  can  be  expressed  in  terms  of  its  phase  advance  0,  related  to  the  reference 
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Figure  /.  Pattern  formation  and  encoding  in  the  olfactory  system  model:  (a)  synthetic  input  pattern  where 
i/O)  are  the  input  signals.  </>,  -phase  shift  of  the  input  signals,  </>w  -  reference  phase  as  in  Equalion  (I);  <b) 
network  of  integrate  and  fire  neurons:  (c)  output  pattern:  u<  -  membrane  potential  of  the  output  neurons. 


phase  (which  corresponds  to  the  zero  time  advance)  as:  <pj  —  4>j/o,  where 
w  —  2n/ T.  In  order  to  distinguish  periodical  time  advances  of  the  input  spikes  from 
the  constant  time  delays  stored  in  the  network,  further  in  the  paper  we  will  call  the 
time  advances  <pj  as  the  phase  shifts  <pj7  related  to  the  reference  phase  <pw,  Thus,  the 
phase  shifts  <pi  of  the  signal’s  spikes  encode  concentrations  of  the  n  corresponding 
constituent  molecules  (2).  An  example  of  input  pattern  is  shown  at  Figure  1(a). 

During  each  cycle  the  spikes  arrive  to  the  neural  layer  with  the  delays  equal  to  their 
phase  shifts  tpj.  Then  the  spikes  acquire  additional  time  delays  d,y  stored  in  the 
synaptic  connections.  So,  the  total  time  delays  of  the  signals  that  arrive  to  tth  neuron 
are  equal  to  ipj  +  d,y. 

Each  neuron  in  the  layer  (Figure  1(b))  is  characterized  by  its  state-membrane 

potential  (i  =  1 . m).  Every  time  a  neuron  receives  a  spike,  its  potential 

Uj  is  increased  by  the  weighted  value  of  that  input  spike:  tv^(r  —  d,y).  At  the  same 
time  the  potential  u,-  is  constantly  decreasing  with  decay  coefficient  k  as  follows: 

+  Y  -  du)  (4) 

/“I 

When  the  potential  of  a  neuron  reaches  its  threshold  value  neuron  fires 

(Figure  1(c)),  and  its  potential  u,-  is  instantly  reset  to  0  as  shown  below. 

Wi(/“)  =  Hthrtth  =>  U/(f+)  =  0  (5) 

The  coefficient  in  (3),  that  is  equal  to  Wthresh/n,  scales  the  value  of  input  spikes.  So, 
only  all  spikes  added  together  are  able  to  increase  the  value  of  potential  to  uthmh-  As 
well  as  the  membrane  potential  ut  is  constantly  decreasing  (4),  the  spikes  that  arrive 
with  big  time  intervals  one  after  another  fail  to  significantly  increase  the  potential's 
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value.  They  have  to  arrive  within  a  narrow  time  interval  (almost  simultaneously)  to 
the  neuron  in  order  to  increase  its  value  to  Uthmh-  Odor  patterns  are  stored  in  the 
network  memory  by  the  delays  d,y  between  inputs  and  the  neurons.  The  delays  are  set 
in  such  a  way  to  make  each  of  n  spikes  arrive  simultaneously  to  the  neuron  if  the 
applied  pattern  is  equal  to  the  pattern  stored: 

di i  =  <P%  ~  tmnUO  (6) 

where  tp §  is  a  phase  vector  of  ith  stored  pattern. 

When  a  stored  odor  pattern  is  applied,  the  total  input  w,yjy(r  -  d,y)  to  the 
neuron  i  has  to  be  equal  to  or  more  than  u,hresh  in  order  to  increase  the  membrane 
potential  ut  by  this  value  and  make  it  fire  (Figure  Ic).  So  the  weights  have  to  be 
determined  as  follows: 


Wij  5* 


uthirsh 

EL.tf'-d*) 


(7) 


The  neuron  fires  simultaneously  with  the  last  arriving  spike,  so  the  phase  shift  of 
the  output  spike  is  equal  to  the  minimal  input  phase  shift  (or  phase  shift  of  the 
weakest  component),  as  shown  in  Figure  2. 

If  the  pattern  applied  is  not  dose  enough  to  any  of  stored  odor  patterns  (of  any 
concentration)  spikes  arriving  with  significant  time  intervals  will  not  make  output 
neuron  fire  because  of  its  exponential  decay.  Output  potentials  in  Figure  1(c)  show  a 
superthreshoid  firing  of  neuron  #5  with  the  phase  shift  and  subthreshold  activity 

of  all  other  neurons.  This  indicates  that  odor  #5  with  some  concentration  is  recognized. 

Due  to  logarithmic  scaling  in  Equation  (2)  the  global  phase  of  the  entire  pattern  is 
shifted  when  the  odor  concentration  is  changed,  while  the  relative  phase  shifts  remain 


(*>) 


Figure  2.  Recognition  of  the  pallems  of  the  same  odor  with  different  concentrations. 
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the  same.  This  makes  pattern  recognition  of  the  model  invariant  to  the  different 
concentrations  of  the  same  odor.  Figure  2  shows  an  example  of  two  patterns  of  the 
same  odor  with  different  concentrations.  Phase  shift  of  the  entire  pattern  (a)  is 
greater  than  of  the  pattern  (b),  though  the  pattern  itself  is  the  same.  Thus,  the 
output  spike  phase  shift  at  (a)  is  greater  than  at  (b).  So,  the  odor  concentration 
of  the  pattern  (a)  is  greater. 

Concentration  of  entire  odor  coyAi  is  defined  as  relative  concentration  of  odor 
applied  to  the  corresponding  odor  stored.  This  concentration  is  decoded  from 
the  output  spike  phase  shift  and  the  minimum  phase  shift  of  the  corresponding 
stored  pattern  min{<p^}  (Figure  2)  by  the  inverse  function  of  Hopfield's  encoding. 

cout,  =  6  exp[(<pom ;  -  min{<p^J)/a]  (8) 


3.  Simulation  Results 

In  our  computational  simulation  the  stimuli  are  4-dimensional  (n  —  4)  concen¬ 
tration  vectors  C  (as  defined  above).  Four  input  neurons  produce  the  spatio-tem¬ 
poral  patterns  that  correspond  to  the  applied  stimuli.  10  output  neurons  (m  =  10) 
correspond  to  10  stored  odors.  The  concentration  of  their  components  varied  from  I 
(threshold  concentration)  to  10.  Concentration  of  the  components  of  the  odors  tested 
varied  from  I  to  70, 

Two  basic  parameters  were  to  be  optimized  during  the  simulation:  exponential 
decay  coefficient  k  (4)  and  corresponding  weights  (7),  The  system  is  quite  sensitive 
to  both  of  them.  With  determined  by  Equation  (7)  with  the  equality  sign,  or  with 
decay  coefficient  k  ;$>  1,  the  system  can  recognize  undistorted  stored  patterns  only. 
When  weights  increase  or  k  decreases  the  network  becomes  more  flexible  (the  recep¬ 
tive  regions  of  the  neurons  get  larger).  However,  the  greater  the  receptive  region,  the 
worse  the  system  accuracy  is.  So,  certain  compromise  had  to  be  found.  For  our 
system  the  highest  success  rate  of  recognition  (shown  in  Table  I)  was  achieved  with 
k  =  6.3,  and  w #  =  1.32,  Vit/  Period  T and  threshold  u^mh  were  arbitrary  selected  as 
follows:  T  =  50,  thresh  —  1*  The  coefficients  at  and  3  define  the  scale  of  the  trans¬ 
formation  (2),  They  were  chosen  as:  oc  =  10,  S  —  lT  that  makes  the  maximum  phase 
difference  of  the  spikes  (that  is  T  —  50)  correspond  to  the  relative  concentration  of 
the  odor  components  equal  to  150. 

Example  of  a  simulated  input  pattern  is  shown  at  Figure  3(a).  This  signal  is  a 
repeating  sequence  of  four  spikes  corresponding  to  applied  odor  {63,96.48,  10} 


Table  /.  Characteristics  of  the  mode!  performance 


Patterns  applied 

5000 

Odors  recognized 

10,3% 

Successful  recognition 

81% 

Incorrect  recognition 

19% 
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Figure  3 .  Simulated  inpul  and  output  patterns:  {a)  input  pattern  of  an  example  odor  {63,96,48.  10);  (b) 
Output  patterns  of  the  1st.  and  Sth  neurons  (both  are  shown  on  the  same  picture).  There  is  a  super- 
threshold  firing  of  the  neuron  #5,  and  subs  thresh  old  activity  of  neuron  Odor  #5  is  recognized. 


which  was  chosen  for  illustration  purpose.  Figure  3(a)  shows  the  activity  patterns  of 
two  neurons.  Only  neuron  #5  reaches  the  threshold  and  fires.  This  indicates  suc¬ 
cessful  recognition  of  odor  #5.  Neuron  #1  does  not  reach  the  firing  level,  so  the  odor 
#]  is  not  recognized. 

In  our  experiment  5000  randomly  generated  test  odors  have  been  presented  to  the 
system.  A  total  of  89.7%  of  them  caused  the  output  neurons  to  produce  the  sub- 
threshold  activity  only  (odors  were  not  recognized).  A  total  of  10.3%  of  the  test- 
stimuli  provoked  superthrcshold  firing  of  output  neurons  (odors  were  recognized).  A 
total  of  81.0%  of  them  were  recognized  correctly  (both  odor  and  concentration), 
with  allowed  concentration  error  equal  to  20%. 

For  the  sake  of  testing  the  proposed  model,  each  applied  odor  is  projected  onto  the 
distance-concentration  coordinate  system.  The  purpose  of  the  test  is  to  determine  if 
our  model  correctly  classifies  an  input  odor  as  resembling  one  of  the  stored  odors,  or 
leaves  the  input  odor  unclassified  if  it  does  not  resemble  any  of  them.  The  resem¬ 
blance  is  expressed  using  a  metric  introduced  in  the  Appendix,  The  stored  odors  have 
a  certain  distribution  in  the  concentration  space.  The  distance  D  between  the  arbi¬ 
trary  input  odor  and  its  closest  neighbor  amongst  the  stored  odors  allows  for  testing 
our  model. 

The  inputs  that  caused  an  output  to  fire  are  indicated  as  points  in  Figure  4(a). 
Majority  of  these  points  are  gathered  in  the  region  of  small  D  with  various  concen¬ 
trations,  This  indicates  that  the  odors  which  caused  an  output  to  fire  were  indeed  close 
to  one  of  the  stored  odor  patterns.  The  concentration  level  of  the  input  odor  does  not 
affect  the  models  ability  to  classify  the  input  as  one  of  the  known  odor  patterns. 

On  the  other  hand,  all  other  input  odors  which  did  not  cause  an  output  to  fire  are 
shown  as  points  in  Figure  4(b).  These  points  tend  to  have  large  values  of  distance 
function  D.  Any  input  odor  that  is  too  far  from  all  of  the  stored  patterns  fails  to 
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Figurt  4.  Distribution  of  the  simulation  results,  (a)  Patterns  correctly  recognized,  (b)  Patterns  that  are  not 
recognized. 


activate  the  model  outputs.  Even  strongly  concentrated  odors  stay  unclassified  if  they 
lack  resemblance  to  one  of  the  stored  patterns* 

In  order  to  compare  Figure  4(a)  and  (b)  a  vertical  division  line  is  drawn.  The  line 
separates  the  classified  points  from  the  unclassified  ones  in  terms  of  the  distance 
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value  D.  Points  to  the  left  from  the  line  manifest  a  resemblance  to  one  of  the  stored 
pattern  odors  and  hence  activate  the  output  of  the  model.  Most  of  the  points  stay  on 
the  right  side  of  the  division  line,  as  no  stored  pattern  claims  such  resemblance. 
Remarkably,  the  division  line  is  vertical,  which  means  that  Hopfieid’s  principle  works 
well  separating  the  odor  'flavors’  from  their  concentrations. 


4,  Discussion  and  Conclusion 

Recognition  of  a  single  pattern  is  only  one  of  the  first  stages  of  sensory  information 
processing.  The  higher-level  problems  of  multi- pattern  processing,  such  as  object 
recognition  in  real  world,  feature  binding  and  segmentation,  object-background 
separation,  attention  focusing,  etc.  are  much  more  difficult  to  model  and  they 
are  far  from  being  handled  by  modem  computational  methods.  However,  all  these 
tasks  are  efficiently  performed  by  cortical  neural  networks  of  animals  and  humans, 
where  different  temporal  correlation  types  are  believed  to  be  the  underlying  prin¬ 
ciple  of  these  abilities  (Campbell  and  Wang,  199S;  Malsburg  and  Buhmann,  1992; 
Malsburg  and  Schneider,  1986). 

Temporal  correlation  plays  an  essential  role  in  olfactory  systems  as  well.  Experi¬ 
mental  results  prove  that  several  odors  in  a  mixture  are  separated  temporally  from 
each  other  at  some  of  the  higher  levels  (Jinks  and  Laing,  1999).  One  of  the  possible 
ways  to  do  such  a  temporal  segregation  is  using  temporal  correlation  and  competi¬ 
tion  of  output  neurons  (or  neural  ensembles)  and  inhibitory  top-down  feedback  to 
input  level  in  order  to  temporally  segregate  recognition  of  different  odors,  suppress 
noise  or  irrelevant  inputs  and  focus  attention  on  the  necessary  odor  (Campbell  and 
Wang,  1998;  Malsburg  and  Buhmann,  1992;  Malsburg  and  Schneider,  1986).  Phase 
encoding,  that  is  a  specific  example  of  temporal  correlation  is  the  basic  principle  of 
our  model  and  we  believe  it  provides  the  base  for  the  solution  of  the  higher  level 
processing  tasks  presented  above. 


Appendix 

This  section  introduces  the  transformation  of  odor  vectors  to  the  distance-concen¬ 
tration  space.  Each  component  cj  of  applied  odor  vector,  C,  is  logarithmically  trans¬ 
formed  as  in  (2)  to  the  corresponding  components  of  phase  vector  q>.  In  the  same 
way  m  stored  vectors  i  =  are  transformed  to  m  vectors  where 

•Pa  =  [<P%J=  I . «}• 

Vectors  <p^  and  <p  are  normalized  in  the  phase  space  in  the  following  way: 
•Pa*  =  \<Paj  ~  ml" \<P%b  j=  l,....  n)  (9) 

<p‘  =  (^-min(^),j=  1 . *} 


(10) 
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Each  applied  and  then  normalized  phase  vector  <p*  =  {<pr  j  =  is  charac¬ 

terized  by  its  distance  D  to  the  closest  of  the  This  distance  D  is  defined  as: 

D  =  min( ||p*  -  (11) 

Thus  D  defines  distance  in  the  phase  space  from  the  vector  applied  to  the  closest  of 
the  stored  vectors* 
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Abstract — The  olfactory  system  is  a  very  efficient  biological 
setup  capable  of  odor  information  processing  with  neural  signals. 
The  nature  of  neural  signals  restricts  the  information  represen¬ 
tation  to  multidimensional  temporal  sequences  of  spikes.  The 
information  is  contained  in  the  interspike  intervals  within  each 
individual  neural  signal  and  Interspike  intervals  between  multiple 
signals.  A  mechanism  of  interactions  between  random  excitations 
evoked  by  odorants  in  the  olfactory  receptors  of  the  epithelium 
and  deterministic  operation  of  the  olfactory  bulb  is  proposed 
in  this  paper.  Inverse  Fro  ben  ms-Perron  models  of  the  bulb's 
temporal  sequences  are  fitted  to  the  interspike  distributions  of 
temporally  modulated  receptor  signals.  Ultimately,  such  pattern 
matching  results  in  ability  to  recognize  odors  and  offer  a  hypo¬ 
thetic  model  for  signal  processing  occurring  in  the  primary  stage 
of  the  olfactory  system. 

Index  Terms — Frobenius  filter,  interspike  intervals,  inverse 
Frobenius  problem,  Markov  process,  odorant  concentration, 
olfactory  bulb,  shift  map,  temporal  sequence. 


I*  Introduction 

LIVING  organisms  perceive  odors  as  sensations  caused 
by  mixtures  of  odorant  molecules.  Such  molecules  ex¬ 
cite  the  olfactory  receptors  to  respond  with  increased  activity 
which  is  then  passed  on  to  the  olfactory  bulb  for  detection. 
Various  odorant  molecules  excite  different  groups  of  receptors. 
A  superposition  of  these  excitations  constitute  the  odor  as 
detected  by  the  olfactory  bulb  [  1  ].  The  relative  concentrations 
of  individual  components  constitute  the  odor  type,  whereas 
the  absolute  concentrations  determine  the  odor  intensity.  The 
olfactory  bulb  has  the  task  of  transforming  the  input  obtained 
from  the  receptors  into  a  set  of  signals  to  be  interpreted  by  the 
brain.  The  capacity  of  a  simple  discriminator  to  distinguish  the 
target  from  background  odorants  has  been  statistically  analyzed 
in  [2]. 

The  continuous  quantities,  such  as  molecule  concentrations, 
cannot  be  directly  represented  by  the  signals  produced  by  bi¬ 
ological  neurons.  Neurons  produce  spikes  and  only  indirectly 
their  presence  or  absence,  or  time  location  may  cany  continuum 
of  information.  The  nature  of  neural  signals  is  assumed  to  have 
the  following  characteristics. 

1)  There  is  no  significance  of  the  shape  of  individual  spikes. 
They  simply  mark  instances  of  time  when  the  neurons  fire. 
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2)  The  signal  is  a  time  sequence  of  spikes.  Spikes  may  occur 
more  or  less  frequently  which  has  an  effect  on  the  average 
value  of  the  signal. 

3)  Spikes  may  occur  in  a  certain  temporal  pattern.  More  pre¬ 
cisely,  the  inter-spike  intervals  may  follow  a  distinct  and 
repetitive  behavior  This  allows  for  code  division  of  infor¬ 
mation  conveyed  by  a  single  signal, 

4)  Two  or  more  signals  may  exhibit  cross  correlation  which 
typically  results  from  synchronization  between  the  signal 
sources.  If  the  synchronized  signals  assume  a  certain 
spatial  distribution,  a  set  of  such  signals  will  manifest  a 
spatio-temporal  pattern. 

The  neural  signals  of  the  olfactory  bulb  representing  the  infor¬ 
mation  about  odors  and  intensities  are  further  interpreted  by  the 
brain.  The  olfactory  bulb  functions  as  the  first  signal  processing 
stage.  In  all  nonbiological  designs  the  first  stage  is  responsible 
for  the  sensitivity  and  noise  performance  of  the  entire  detection 
system.  The  same  should  hold  true  in  case  of  the  olfactory  bulb, 
A  very  detailed  investigation  of  neuronal  noise  and  spike  prop¬ 
agation  can  be  found  in  [3]. 

The  goal  of  this  article  is  to  identify  the  simplest  method 
of  encoding  odor  information  in  temporal  sequences.  The 
input-output  interactions  between  temporal  sequences  can  lead 
to  an  odor  detection  and  encoding  mechanism  in  the  olfactory 
bulb. 

II-  Temporal  Modulation  1 

The  very  input  of  the  olfactory  system,  the  epithelium,  pro¬ 
duces  an  enormous  number  of  signals.  Receptors  are  hard- wired 
to  detect  specific  odor  components  and  are  uniformly  distributed 
in  the  epithelium.  The  odor  information  is  therefore,  spatially 
distributed  across  the  epithelium  and  is  assumed  to  have  no  tem¬ 
poral  dependency.  Every  odor  and  concentration  can  be  repre¬ 
sented  by  its  "black  and  white  photo”  in  which  the  gray  levels 
of  pixels  encode  spiking  activities  of  the  receptors.  In  this  paper, 
the  odor  information  is  assumed  to  be  spatially  distributed  and 
static,  although  there  is  a  strong  evidence  of  various  significant 
aspects  of  the  inhale-exhale  rhythm  and  the  impulse  response 
of  the  olfactory  bulb  [5]. 

No  temporal  coding  of  information  is  performed  by  the  indi¬ 
vidual  receptors.  Simply,  the  more  molecules  are  present  at  the 
docking  sites  of  the  receptors,  the  more  frequent  their  spiking. 
Based  on  response  measurements  and  fitting  of  concentration- 
response  curves  presented  in  [6]  and  [7],  the  spiking  frequency 
/  of  the  receptors  has  an  asymptotic  dependency  on  the  odorant 
concentration  c  (molarity).  When  the  odorant  concentration  is 

’Term  temporal  modulation  is  adopted  from  [4], 
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at  the  lowest  detectable  level  c  =  c,,  the  receptor  fires  at  the 
very  slow  rate  of  spontaneous  activity.  When  the  concentration 
grows  infinitely  large,  the  frequency  reaches  saturation  at  the 
maximum  firing  rate  of  fm.  The  curve  f(c)  fits  the  following 
definition: 

2 

f  =  “/m  axctan  0(c  -  cf).  (1) 

The  slope  factor  0  is  expressed  in  terms  of  the  dynamic  range  Ac 
defined  as  the  odor  concentration  at  which  the  frequency  reaches 
80%  maximum,  /{c*  + Ac)  =  0 ,8/m.  Given  the  dynamic  range, 
the  slope  factor  can  be  determined  as  0  =  tan(0.87r/2)/Ae. 

Concentration  c,  used  in  [6]  and  [7],  is  a  logarithmic  quantity 
related  to  the  odorant  molarity  c  =  logm,  with  m  in  mol/L 
The  investigated  odorants  were  anisole  (ANI),  camphor  (CAM), 
isoamyle  acetate  (ISO),  and  limonene  (LIM).  The  curve  fitting 
resulted  in  the  following  parameters  for  each  odorant  [6]: 


ANI 

CAM 

ISO 

LIM 

fm 

11  Hz 

15  Hz 

11  Hz 

8  Hz 

Ct 

-6.7 

-8.6 

-7.0 

-7.7 

l 

1.1 

1.1 

0.5 

0.3. 

Remarkably,  linear  curves  are  obtained  if  instead  of  spiking 
frequency  /,  the  interspike  intervals  r  -  l/f  are  graphed 
versus  the  reciprocal  of  molarity,  referred  to  as  sparsity  s .  Since 
different  odorants  may  have  largely  different  molarity  threshold 
levels  a  =  logmt,  the  reciprocal  of  the  incremental  molarity 
$  —  1  /(m  —  m()  rather  than  the  absolute  value  can  be  used 
in  the  joint  graph  for  various  odorants.  The  parametric  repre¬ 
sentation  of  relationship  (1)  in  the  new  coordinates  (£,r)  for 
the  introduced  odorants  is  shown  in  Fig.  I.  The  horizontal  and 
vertical  axes  are  the  incremental  sparsity  s  and  the  interspike 
interval  r  expressed  in  terms  of  molarity  m  as  follows: 


r(m)  =  arctan  0  log  .  (3) 

As  can  be  seen  in  the  figure,  diluting  the  odorant  in  the  air 
increases  the  imerspike  intervals  at  an  approximately  constant 
rate.  This  may  be  regarded  as  temporal  modulation  with  the  con¬ 
version  gain  G  =  dr/ds%  which  is  the  slope  of  the  line.  The  left 
side  of  each  curve  corresponds  to  the  receptor  saturation  region. 
By  extrapolating  the  curves  to  the  intersections  with  the  ver¬ 
tical  axis,  a  minimum  interval  tq  for  each  receptor  type  can  be 
found  of  value  roughly  around  i  00  ms.  This  minimum  interval 
may  be  regarded  the  refractory  period  of  the  receptor  With  just 
two  parameters  r0  and  G  for  each  receptor  type,  the  temporal 
modulation  illustrated  in  Fig.  I  can  be  readily  described  using 
first-order  approximation: 

r  =  r0  -f  Gs.  (4) 

The  approximation  can  be  validated  only  within  the  dy- 
namic  range  of  the  receptor,  that  is,  outside  the  saturation 
5  >  l/(mt10Ac). 


Ill  Odor  Characterization  With 
Interspike  Distributions 

An  odor  is  a  superposition  of  a  number  of  basic  odorants. 
The  concentration  information  is  temporally  modulated  at  the 
glomerular  inputs  of  the  olfactory  bulb,  therefore,  the  percep¬ 
tion  of  odor  intensity  must  be  related  to  the  imerspike  intervals. 
Increasing  the  odor  intensity  shortens  the  intervals  at  different 
rates  for  each  basic  odorant  due  to  the  differences  in  their  con¬ 
version  gains.  This  provides  some  explanation  why  responses  of 
the  mitral  outputs  can  be  different  for  the  same  odor  at  different 
intensities  [8], 

In  the  glomerular  layer  the  enormous  number  of  inputs 
converge  into  much  less  dimensional  connections  to  the  mitral 
cells.  The  glomeruli  are  also  highly  interconnected  between 
each  other  via  periglomerular  imemeurons  [9].  Both  inhibitory 
and  excitatory  connections  are  present  within  the  glomeruli 
which  indicates  that  a  winner-take-all  mechanism  could  be  in¬ 
volved  before  the  input  to  the  mitral  cells.  The  presence  of  such 
a  mechanism  would  enable  arranging  of  the  input  interspike 
intervals  into  distributions  statistically  representing  the  odor 
information. 

Let  R  be  the  number  of  ail  types  of  receptors  in  the  epithe¬ 
lium.  This  makes  R  also  the  number  of  distinct  basic  odorants, 
the  basis  for  the  odor  space.  Suppose  the  first  four,  out  of  all  R 
odorants,  are  the  ones  shown  in  Fig.  I  .  An  odor  at  a  given  in¬ 
tensity  can  be  uniquely  represented  by  a  vector  s  of  sparsities  of 
the  basic  odorants.  For  instance,  an  odor  created  by  mixing  0.5 
mol  of  CAM  and  0.75  mol  of  LIM  with  100  Ml  of  air  would  be 
represented  by  vector  a  =  (oof  50,  «>,  75,  oo, . . . ,  oo)  Ml/mot. 
Vector  2s  would  represent  the  same  mixture  diluted  in  twice  the 
amount  of  air.  In  general,  an  odor,  as  seen  by  the  epithelium, 
is  s  -  ($its2l  Terms  “vector11  and  “basis"  are  under¬ 

stood  to  be  suitable  ways  to  arrange  numbers  rather  than  the 
strictly  defined  terms  used  in  linear  spaces, 

A  much  more  compressed  way  to  describe  odors  is  through 
distributions  of  interspike  interval  probabilities.  This  formalism 
may  also  be  more  relevant  to  the  signals  presented  to  the  mitral 
inputs.  Let  the  interspike  intervals  be  quantized  into  TV  ranges 
with  cutoff  rmax.  Interval  rmiLX  is  considered  to  be  a  borderline 
between  evoked  and  spontaneous  activity  of  the  receptors,  A 
single  neural  signal  can  represent  an  odor  with  the  interspike 
interval  r  probability  distribution  p  =  *  *  ,pjv)t  which 

is  formally  a  vector  of  probabilities 

p  =  <  r  -  7?Tmax)  s  if  n  <  TV 

lPr{rmax<r),  ifn  =  TV. 

The  quantized  representation  of  the  interspike  interval  distri¬ 
bution  is  chosen  because  it  is  more  suitable  for  numeric  com¬ 
putations  than  the  probability  density  function.  A  satisfactory 
approximation  of  continuum  can  be  achieved  provided  that  TV 
is  large  enough. 

Suppose  the  0,5  mol  of  CAM  and  0.75  mol  of  LIM  mixture 
with  100  Ml  of  air,  indicated  by  the  filled  squares  in  Fig.  I,  is 
presented  to  the  olfactory  epithelium.  Two  kinds  of  receptors 
would  be  activated,  each  responding  with  spikes  separated  by 
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lOOMI/mol  200 Ml/mol  300MI/mol 


Fig.  2.  Odor  composed  of  0.5  mol  of  CAM  odoram  and  0.75  mol  of  LIM 
odoram  mixed  with  100  Ml  of  air  (filled  bars)  and  then  diluted  in  additional 
100  Ml  of  air  (empty  bars).  The  bars  have  width  rmi^/.V  and  represent 
probabilities  of  respective  time  intervals  as  defined  by  (5).  Example  with 
N  =  20  and  Tmmx  =  350  ms  shown. 


^■8-  y  Interspike  intervals  r  of  receptor  firing  versus  incremental 

sparsity  a  of  odorant.  Six  points  on  each  curve  correspond  to  the  following 
logarithmic  concentration  levels  (lefHo-righi);  1.2 Ac,  Ac,  0.8 Ac,  0.6A c, 
0.4 Ac,  and  0.2 Ac  above  c* .  The  axes  units  are  ms  and  Ml/mol, 


invariant  distribution  equal  p*  could  serve  as  a  first-order  ap¬ 
proximation  of  a  dynamical  system  for  that  odor.  Let  N  x  N 


matrix  P  be  the  transition  matrix  of  the  Markov  process 


roughly  90-ms  intervals  and  1 75-ms  intervals,  respectively,  ac¬ 
cording  to  Fig.  L  Suppose  also  that  there  is  twice  as  many  LIM 
receptors  as  those  detecting  CAM.  In  the  winner-take- all  com¬ 
petitions,  the  LIM  receptors  would  have  a  better  chance  passing 
its  signal  compared  to  the  CAM  receptors.  The  described  odor 
is  represented  by  the  filled  bars  of  interspike  interval  probabil¬ 
ities  in  Fig.  2,  The  probability  of  the  175  ms  LIM  intervals  is 
twice  the  probability  of  the  90  ms  CAM  intervals. 

Suppose  further  that  the  same  odor  mixture  is  now  diluted 
in  twice  the  amount  of  air.  This  doubles  the  sparsity  of  the 
odor,  hence,  increasing  the  imerspike  intervals  of  both  odor¬ 
ants  present  in  the  mixture.  The  diluted  odors  are  represented  in 
Fig.  2  by  the  empty  bars  of  probabilities.  Now  the  LIM  intervals 
are  about  220  ms  and  the  CAM  intervals  increased  to  about  100 
ms  without  a  change  to  the  probability  levels.  Note  that  the  two 
odorants  have  different  conversion  gains  and  modulate  the  tem¬ 
poral  intervals  at  different  rates.  As  the  odor  intensity  changes* 
this  changes  the  probability  pattern,  A  different  hypothesis  of 
time  advance  modulation  where  the  resulting  partem  is  invariant 
under  the  concentration  level  was  introduced  in  [10]  and  leads 
to  functional  models  [11],  However,  the  neurophysiological  ev¬ 
idence  suggests  that  the  patterns  of  bulbar  activity  actually  do 
change  when  varying  the  odor  intensity  [8], 

The  signal  processing  occurring  between  the  mitral  cells  and 
glomerular  layer  is  a  dynamical  process.  The  information  is  em¬ 
bedded  in  the  time  realizations  of  signals.  It  may  be  retrieved 
only  through  observation  of  these  signals  for  a  period  of  time. 
The  probability  distribution  of  the  interspike  intervals  may  be 
retrieved  by  statistically  analyzing  the  neural  signals.  Likewise, 
a  simple  stochastic  process  can  be  modeled  to  have  the  statis¬ 
tical  properties  representing  a  given  odor  through  the  probability 
distribution, 

Suppose,  in  steady-state  after  all  the  transient  response  has 
vanished  the  odor  is  represented  by  the  probability  distribution 
p*  defined  according  to  (5),  A  Markov  process  [12]  with  the 


p(k  +  l)  =  Pp(k)>  (6) 

Also,  let  the  process  converge  to  p*  in  a  sense  that 
P m  =  lim^^p^)  for  almost  every  initial  distribution 
p(Q).  The  invariant  distribution  is  the  eigenvector  of  transition 
matrix  P  associated  with  the  unit  eigenvalue;  p*  =  Pp V  In 
this  respect,  the  Markov  process  is  a  dynamical  system  in 
probabilistic  space  S  =  {p  e  (0;  l]^  |  pn  -  1}  with  a 
stable  fixed  point  p*.  Further  on,  space  $  will  be  referred  to  as 
the  odor  space. 

Consequently,  an  odor  may  be  associated  with  an  operator 
P  :  S  S  in  the  odor  space.  The  odor  itself  is  the  stable  fixed 
point  of  the  operator.  There  is  a  benefi  t  of  such  a  representation 
of  odors.  Operator  P  defines  an  odor  indirectly  through  a  defi¬ 
nition  of  a  dynamical  system.  It  is  easy  and  natural  to  generate 
realizations  of  neural  signals  using  such  operators,  which  is  suit¬ 
able  in  the  modeling  effort.  There  are  many  operators  that  have 
the  same  invariant  distribution.  Hence,  the  same  odor  informa¬ 
tion  may  be  redundantly  embedded  in  many  different  processes. 
Another  approach  to  representing  odors  with  Markov  processes 
is  presented  in  [13], 

Formally,  a  realization  of  the  introduced  Markov  process  is  a 
sequence  of  imerspike  intervals  {r*  }.  Define  the  interval  range 
to  be  Tn  —  [(n  —  I/AT)rmax;  Tmax)  if  n  <  N>  and  T/v  — 

;  oo)  otherwise,  where  the  interval  range  index  n  is  defined 
in  the  same  manner  as  in  (5).  For  the  sake  of  modeling  through, 
optimization,  a  particular  operator  P  may  be  developed  to  have 
p*  as  its  invariant  distribution  of  interspike  intervals  over  time. 
Denote  the  elements  of  the  operator  by  f  so  that  P  —  [F^], 
where  i  and  j  are  the  row  and  column  indexes.  Number  is 
the  probability  that  in  the  Markov  process  (6)  an  interval  from 
the  range  T*  will  follow  the  interval  from  the  range  7}: 

p  _  Fr(rfc+1  €  Tj  and  rk  E  Tj) 
lJ  Pr(T*  e  Tj) 


(7) 
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There  is  no  dosed- form  formula  for  selecting  Pi}s  for  a  given 
p  .  However,  starting  with  some  random  PtjS,  an  optimization 
algorithm  can  be  used  to  find  the  PijS  as  the  minimum  of  a 
suitable  cost  function.  Since  ail  P^s  are  probabilities,  they  must 
be  numbers  in  the  unit  segment  from  0  to  1 .  This  fact  allows  for 
constructing  one  of  the  components  to  be  included  in  the  cost 
function,  namely,  the  unit  segment  potential.  For  each  number 
Pij  > a  potential  function  v(P,j ),  shown  in  Fig.  3,  describes  how 
distant  P^  is  from  the  unit  segment 


v(Pij) 


(2 Pij  - 1)2 
l  +  (2Pij-  !)-«■ 


(8) 


Function  v  attains  the  minimum  in  the  middle  of  the  unit  seg¬ 
ment  and  is  maximally  flat  within  the  segment.  The  maximally 
flat  approximation  [14]  with  a  rational  function  is  chosen  to  fa¬ 
cilitate  the  optimization  process.  The  partial  costs  v(Pij )  sum  up 
to  the  cost  function  component  EV(P)  responsible  for  keeping 
all  the  entries  of  P  within  the  unit  segment 


w=EEw  (9) 

*=i  j=i 

Operator  P  is  a  probabilistic  matrix  in  a  sense  that  all  its 
column  vectors  are  normalized  probability  distributions.  There¬ 
fore,  the  column  sums  of  P  must  sum  up  to  L  Another  cost 
function  component  EC(P)  measures  the  deviation  from  this 
requirement 


n  /  n  \  2 

E'(P)  -  £  1  -  .  (10) 

i=i  V  t=i  / 

If  operator  P  is  a  well  defined  transition  matrix  of  Markov 
process  (6),  then  the  cost  sum  Ev  (P)+ J5C(P)  is  low  and  close  to 
its  minimum  attainable  value.  The  goal  of  the  cost  minimization 
procedure  is  to  develop  operator  P  with  the  constraint  that  p*  is 
it's  principal  eigenvector  associated  with  eigenvalue  1 .  To  sim¬ 
plify  the  operator  synthesis,  matrix  P  will  be  assumed  to  be  di- 
agonalizable:  P  =  BAB-1.  The  diagonal  matrix  A  =  diag(A) 
is  composed  of  N  eigenvalues  A  -  (Ai  ,  A2,  * .. ,  XN)  of  P. 
Let  Aj  =  L  The  convergence  rate  of  the  dynamical  system  (6) 
heavily  depends  on  the  radius  of  the  remaining  A,s  for  i  >  1. 
Operator  P  is  synthesized  with  random  A*s,  for  i  >  1,  with 
the  assumption  that  |A(|  <  r  <  1  and  the  radius  r  is  kept  low 
to  improve  the  convergence  rate.  In  the  numerical  experiment  r 
was  selected  to  be  equal  to  0,2.  Operator  P  is  diagonal  in  the 
basis  constructed  with  the  column  vectors  of  B,  Since  Ai  =  1, 
the  first  column  vector  of  B  is  p* .  More  precisely,  =  p* 
for  j  =  1.  All  other  entries  Bij ,  for  j  >  1,  are  variables  in  the 
optimization  process.  Their  initial  values  are  selected  randomly 
from  the  uniform  distribution  in  the  range  (-1;  1).  Final  matrix 
B  is  found  using  an  optimization  algorithm  to  minimize  the  cost 
function  as  in  the  following  expression: 


min  [EV(P)  +  EC(P)1. 

lj>l 


(11) 


p 

U 

Fig.  3.  Unit  segment  potential  function.  In  the  total  cost  function,  numbers 
Pii  €  [0;  I)  contribute  much  less  than  numbers  PtJ  outside  this  range 
Minimisation  of  w  will  attract  all  P^  s  toward  the  inside  of  the  unit  segment 

zeroing  of  negative  values  and  normalization  of  columns.  The 
principal  eigenvector  of  P  is  not  very  sensitive  to  such  trimming 
of  P. 

IV.  Embedding  Distributions  in  Temporal  Sequences 

As  illustrated  in  the  example  shown  in  Fig.  2,  the  probabilistic 
representation  of  odors  and  intensities  fits  well  the  random  na¬ 
ture  of  excitations  received  from  the  olfactory  epithelium.  The 
Markov  model  is  also  a  natural  candidate  for  a  simple  approx  i- 
mation  of  the  dynamics  behind  spike  interactions  driven  by  the 
receptors.  The  olfactory  bulb,  however,  should  be  considered  a 
deterministic  system  which  has  no  random  variables  other  than 
the  input  received  from  the  epithelium.  Moreover,  the  olfac¬ 
tory  bulb  is  capable  of  self-excitatory  activity  in  response  to  the 
input.  This  may  be  the  factor  contributing  to  both  high  sensi¬ 
tivity  and  high  selectivity  of  the  sense  of  smell  [15].  From  this 
perspective,  it  seems  reasonable  to  regard  the  olfactory  bulb  as 
an  active  medium  rather  than  a  passive  relay  of  receptor  signals. 
The  olfactory  bulb  actively  produces  firing  activity  in  response 
to  the  receptor  signals  [16], 

A  sequence  of  interspike  intervals  complying  with  a  given 
interval  distribution  can  be  generated  in  a  deterministic  dynam¬ 
ical  system.  The  simplest  such  system  is  a  one-dimensional  map 
constructed  by  solving  the  inverse  Frobenius-Perron  problem 
[  1 7].  The  overall  goal  of  the  search  for  a  sequence  generator  is 
to  be  able  to  represent  the  odor  information  by  a  distribution  of 
interspike  intervals.  A  simple  shift  map  can  be  constructed  di¬ 
rectly  from  the  probabilistic  operator  used  in  approximation  (6), 
as  described  in  detail  in  [18]. 

First,  a  piecewise  linear  map  /  :  [0;  N ]  ->  [0;  N],  [0;  N]  c  R 
is  derived  from  probabilities  included  in  the  operator  P 


ifi  ~Yl  Pmi<X<j  +  Pij  -  Prny  (12) 

m=l  m®I 


V(P  } 

U 


The  minimized  solution  oftentimes  needs  a  final  touch  to  make 
sure  that  P  has  no  negative  entries,  no  entries  greater  than  one 
and  that  column  sums  of  P  are  indeed  ] .  This  can  be  done  by 


As  shown  in  Fig.  4,  map  /  is  composed  of  N 2  linear  segments 
corresponding  to  N2  numbers  Pijt  The  slopes  of  the  segments 
are  simply  If  Pij.  To  evaluate  f(x ),  the  pair  of  indexes  i  and  j 
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appropriate  for  a  given  x  needs  to  be  identified  using  the  condi¬ 
tion  provided  by  (12). 

By  scaling  the  domain  and  range  of  map  /.  the  dynamical 
system  generating  temporal  sequence  {rfc}  can  be  defined  with 
the  help  of  shift  map  h 


k(r)  = 


N 


1  \  Tmmc  / 


(13) 


Regardless  of  the  initial  condition  chosen,  subsequent  map¬ 
pings  with  h  will  determine  sequence  {rk}  whose  distribution 
of  values  converges  to  the  invariant  distribution  of  process  (6)* 
A  deterministic  dynamical  system 


n+i  =  A(t*)  (14) 

may  be  regarded  as  a  generator  of  realizations  of  neural  signals 
for  a  given  distribution  of  interspike  intervals- 
A  numeric  example  of  shift-map  synthesis  is  shown  in 
Fig.  5.  Three  interspike -interval  distributions  p*AJ  p*B ,  and  p*c 
are  selected  randomly  to  characterize  three  hypothetical  odors  A, 
Bt  and  C.  The  bars  represent  probabilities  of  respective  time 
intervals  as  defined  by  (5).  The  horizontal  axis  is  normalized 
such  that  interval  NrTT1AX/(N  —  1)  corresponds  to  1.  Shift 
maps  hAr  hB ,  and  hc  are  evaluated  for  the  example  odors 
and  shown  in  the  middle  row  of  figures  also  in  time-normalized 
coordinates.  The  maps  have  N  disconnected  branches,  however, 
vertical  lines  connecting  the  branches  are  added  to  enhance 
the  graphs.  Starting  with  a  random  initial  interval,  each  map 
iterated  3000  times  according  to  (14)  produced  a  temporal 
sequence.  The  sequences  are  shown  in  the  bottom  row  of 
figures.  Each  interval  in  a  sequence  is  indicated  as  a  point 
whose  vertical  coordinate  is  the  normalized  time.  This  way 
the  density  of  points  reflects  the  original  distribution  plots  if 
rotated  clockwise  by  90°. 


V,  Frobenius  Filter  for  Temporal  Sequences 

It  is  broadly  accepted  that  the  olfactory  bulb  provides  sup¬ 
port  for  a  pattern  recognition  mechanism  for  odor  detection  and 
classification.  Not  all  of  the  recognition  is  taking  place  there, 
but  definitely  the  process  is  initiated  in  the  olfactory  bulb.  As¬ 
suming  that  the  temporal  sequences  of  interspike  intervals  are 
earners  of  odor  information,  an  implementation  of  signal  pro¬ 
cessing  system  (14)  can  be  proposed.  Ultimately,  the  goal  is 
to  demonstrate  usefulness  of  the  proposed  mechanism  in  odor 
recognition. 

The  signal  processing  scheme  shown  in  Fig.  6  will  be  refen-ed 
to  as  the  Frobenius  filter.  The  input  to  the  filter  is  a  temporal  se¬ 
quence  whose  interspace  intervals  are  determined  by  the  random 
variable  rin  with  values  governed  by  the  probability  distribution 
Pin  defined  as  in  (5),  Distribution  pin  characterizes  an  odor. 

The  Frobenius  filter  is  simply  a  shift  map  with  the  feedback 
loop  controlled  by  a  random  switch.  The  switch  operation  is 
described  by  a  two- valued  stochastic  process  f  :  {0, 1}  x  N  -* 
R.  The  filter  is  producing  time  intervals  based  on  the  switch 
position.  At  every  interval  k%  the  switch  position  depends  on  the 
value  of  &  governed  by  probabilities 

(15) 

(16) 


Hg,4.  Example  of  a  piecewise  linear  shift  map /(i).  Function  /  is  composed 
Of  iV  continuous  branches  fs  \  f>  -  1;  j)  [0;  AT],  If  *  is  chosen  randomly 
from  the  uniform  distribution  over  the  range  (0;  .V).  the  conditional  probability 
€  (*  “  1;  0  given  the  fact  *  6  (j  -  1; j)  can  be  evaluated  by 
-  fl/j  (l1  “  1  i  »])||.  Example  with  N  -  3  shown 


where  c  €  [0;  1]  is  a  constant  parameter.  When  =  1,  the  filter 
is  receiving  the  input  rin.  The  opposite  position  of  the  switch 
(?*  =  0)  the  shift  map  determine  the  output  time  interval 
based  on  the  previous  interval  as  in  (14).  The  overall  filter  equa¬ 
tion  reads 


T-fe+i  =*  h  [&rm  +  (1  -  &)t*] .  (17) 

The  notion  of  the  switch  is  an  attempt  to  model  a  competition 
between  the  input  and  the  feedback.  Its  random  operation  is  in¬ 
herited  from  the  random  nature  of  the  input  temporal  sequence. 

The  three  shift  maps  hA,  hs,  and  he,  introduced  in  Fig.  5, 
are  used  to  illustrate  the  function  of  the  Frobenius  filter.  Each  of 
the  shift  maps  was  stimulated  at  the  input  by  values  generated 
by  probability  distributions  p^,  and  p’c  representing  three 
different  odors.  Fig.  7  shows  all  possible  input-filter  combina¬ 
tions  arranged  in  the  following  nine  pairs: 

(Pa>^a)  (pb,/ia)  (pci  hA) 

(Pa>  ha)  (Pa,  hs)  (Pc-^b) 

(Pa.Ac)  (Pb,M  (Pc.fcc).  (18) 

In  each  instance,  K  =  20000  values  random  values  rin  were 
drawn  from  the  input  distribution  and  applied  with  probability 
c  =  0.5  to  the  filter.  The  values  drawn  were  also  sorted  in  the 
ascending  order  and  stored.  In  the  sorted  input  sequence  {u,} 
the  following  property  holds:  i  <  j  =>u,-  <  Uj.  When  plotted, 
the  graph  of  the  sorted  sequence  would  resemble  the  shape  of 
the  cumulative  distribution  function  of  the  random  variable  rin. 

The  realization  of  the  sequence  {r*}  generated  by  the 
Frobenius  filter  for  K  iterations  were  also  sorted  in  the  same 
manner.  The  sorted  output  sequence  {£;}  was  then  compared 
to  the  sorted  input  sequence  in  Fig.  7.  More  precisely,  the 
graphs  in  the  figure  are  the  sequences  of  quadratic  distances 
{(ui  -  fi)2}  in  each  of  the  nine  instances.  The  horizontal  line 
is  the  mean-square  value  of  the  distance  (u,  -  U)2.  As  seen 
in  the  figure,  the  input-output  sequences  generated  in  pairs 


Pr(£*  =  1)  =  c 
Pr(6t  =  0)  =  1  -  c 
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the  lop  row!  The  corresponding  shift  m^s^dcUsfr^wl^ofwSues  ge Termed  remMrd  Tea  ucncL^-e  'the  ~  20rePrese"tii^ thrtc  different  odors  are  shown  in 
to  «.  range  of  (0;  I ).  Bach  graph  in  the  hot, ore  row  contains  3000  points 


Pt  =  Jc 

Shiftmap 

Pr  =  c 

Fig.  6.  Frobenius  filter  is  a  shift  map  with  input.  Either  the  input  interval  rin  or 
the  present  output  interval  rk  is  transformed  into  the  neat  output  interval  rt+1, 

(Pa>  ^a),  (Ps,  /is),  and  (p'c%  he)  are  synchronized  in  a  sense 
that  the  quadratic  distance  between  input  and  output  interval 
distributions  is  small.  The  distances  in  all  the  other  pairs  are 
significantly  larger.  By  detecting  low  distance  between  the  input 
and  the  output  of  the  filter,  an  odor  recognition  mechanism  can 
be  devised. 

Two  examples  of  realizations  of  the  input  and  output  temporal 
sequences  are  shown  in  Figs.  8  and  9.  The  proposed  mecha¬ 
nism  uses  a  pattern  matching  phenomenon  which  signals  suc¬ 
cessful  detection  as  a  decreased  distance  between  parameters  of 
the  input  and  the  output  neural  signals,  The  pattern  matching 
is  not  based  on  coherence  of  the  two  signals,  As  shown  in  the 
figures,  no  similarity  in  realizations  of  the  input  and  the  output 
can  be  observed  in  either  the  matched  or  unmatched  odor-filter 
pairs.  In  case  of  the  matched  pair,  the  similarity  is  in  the  statis¬ 
tical  properties  of  the  input  and  the  output  signals. 

Value  c  =  0.5  used  in  the  experiment  is  the  midpoint  of  the 
probability  range.  The  value  of  c  may  be  selected  to  optimize 
the  performance  with  respect  to  different  numbers  of  odors  and 


filters.  The  extremes  correspond  to  no-input  (c  =  0)  and  the 
open-loop  (c  =  I)  conditions.  In  this  regard,  c  represents  the 
strength  of  the  input  coupling.  With  no  input,  the  entire  system 
becomes  a  simple  pattern-matching  mechanism.  In  the  open- 
loop  condition,  the  output  is  driven  by  the  input  through  the  filter 
mapping  h.  In  both  cases  if  the  input  is  a  stationary  process, 
the  filter  output  is  also  stationary  and  the  distance  u,  —  l,  can 
be  evaluated.  However,  there  will  be  no  gain  resulting  from  the 
synchronizing  effect  imposed  by  the  switch. 

Vi.  Conclusion 

The  details  of  how  the  cells  of  the  olfactory  bulb  could  en¬ 
code  the  information  in  the  way  described  in  this  paper  are  not 
discussed  here.  The  goal  set  for  this  work  was  to  describe  the 
simplest  method  of  encoding  information  in  temporal  sequences 
and  show  the  input-output  interactions  which  can  lead  to  an  odor 
detection  and  encoding  mechanism.  All  the  computations  are 
very  simple.  No  memory  is  necessary,  only  the  last  time  interval 
is  locally  kept  in  the  evaluation  of  the  next  time  interval  in  the 
output  sequence.  Actual  neurons  are  capable  of  performing  such 
a  storage  with  their  inherent  leaky  integration. 

The  computational  complexity  of  the  model  depends  on  the 
resolution  N  of  the  interspike  interval  distributions.  It  deter¬ 
mines  the  dimensionality  of  the  Markov  transition  matrix.  Mini¬ 
mization  (1 1)  is  the  most  time-demanding  computation  involved 
in  the  approach  and  takes  a  significant  amount  of  time  to  eval¬ 
uate.  The  minimization  may  be  regarded  as  the  learning  process. 
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?f;  •l.,:.N.Ule.rn!itajlCeS1.0f  “  Fn*enius  filler  stimulated  by  an  input  distribution  for  20  000  iterations.  Quadratic  distance  <f(  =  (u.  -  t.P  between  cumulative 
distributions  of  mterspike  intervals  at  input  «,  and  output  f.  of  the  filter  shown.  The  graphs  are  arranged  according  to  (18).  cumulative 
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Fig,  8.  Realizations  of  the  input  {top)  and  output  (bottom)  temporal  sequences 
for  a  matched  pair  { p \ ,  h* )-  A  fragment  containing  100  spikes  shown. 
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hig.  9.  Realizations  of  the  input  (top)  and  output  (bottom)  temporal  sequences 
for  unmatched  pair  (p^,  h*  ).  A  fragment  containing  100  spikes  shown. 

The  shapes  of  the  shift  maps  shown  in  Fig.  5  are  not  crucial 
for  in  the  operation  of  the  Frobenius  filters.  The  proposed  shapes 
are  just  samples  of  infinite  possibilities  chosen  for  mathematical 
simplicity.  In  fact*  any  mapping  that  is  mixing  and  expanding 
[19]  can  be  used  in  the  Frobenius  filter.  Such  mappings  develop 
continuous  invariant  distributions  and  guarantee  ergodicity  of 
the  temporal  sequence  in  a  sense  that  the  invariant  distribution 
is  reachable  from  any  initial  condition.  Temporal  sequences  at 
the  output  of  the  filter  have  a  very  strong  ability  of  encoding 
information  in  the  time  scale.  It  is  sufficient  to  isolate  just  a  few 
consecutive  spikes  to  be  able  to  identify  the  shift  map  which 
generated  the  spikes  and  effectively  identify  the  encoded  odor. 
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Abstract 

The  ability  of  bifurcating  processing  units  and  their  networks  to  rapidly  switch  between  different  dynamic  modes  has  been  used  in  recent 
research  efforts  to  model  new  computational  properties  of  neural  systems.  In  this  spirit,  we  devise  a  bifurcating  neuron  based  on  control  of 
chaos  collapsing  to  a  period-3  orbit  in  the  dynamics  of  a  quadratic  logistic  map  {QLM).  Proposed  QLM  3  neuron  is  constructed  with  the  third 
iterate  of  QLM  and  uses  an  external  input,  which  governs  its  dynamics.  The  input  shifts  the  neuron's  dynamics  from  chaos  to  one  of  the  stable 
fixed  points.  This  way  the  inputs  from  certain  ranges  (clusters)  are  mapped  to  stable  fixed  points,  while  the  rest  of  the  inputs  is  mapped  to 
chaotic  or  periodic  output  dynamics,  ft  has  been  shown  thai  QLM 3  neuron  is  able  to  learn  a  specific  mapping  by  adaptively  adjusting  its 
bifurcation  parameter*  the  idea  of  which  is  based  on  the  principles  of  parametric  control  of  logistic  maps  [Proceedings  of  the  International 
Symposium  on  Nonlinear  Theory  and  its  Applications  (NOLTA*97)t  Honolulu,  HI,  1997;  Proceedings  of  SP1E,  2000]  Learning  algorithm 
for  the  bifurcation  parameter  is  proposed*  which  employs  the  error  gradient  descent  method. 

©  2003  Elsevier  Ltd.  All  rights  reserved. 
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1.  Introduction 

Limitations  of  the  static  nature  of  artificial  neural 
networks  (ANN)  stimulate  investigation  of  biologically 
motivated  neuron  models  with  inherent  dynamics*  such  as 
bifurcating  and  chaotic  neurons  (Farhat*  2000;  Farhat  & 
Eldefrawy*  1992;  Farhai*  Lin*  &  Eldefrawy,  1994;  Holden, 
Hyde,  Muhamad,  &  Zhang*  1992).  In  the  networks 
composed  of  such  neurons  information  is  processed  by 
convergence  not  only  to  a  fixed  point,  but  also  to  a  limit 
cycle  or  chaotic  attractor  (Farhat*  2000;  Farhat,  Lee*  &  Ling* 
1998;  Hirsch  &  Baird.  1995;  Lee  Sl  Farhat*  2001), 

Artificial  neurons  commonly  used  to  mimic  dynamics  of 
biological  neurons  are  simplified  versions  of  the  Hodgkin- 
Huxley  type  model  (HHM)*  such  as,  for  instance*  integrate- 
and-fire  neurons*  These  models  demonstrate  very  rich 
dynamics*  with  a  variety  of  bifurcations  and  chaotic 
phenomena  (Farhai  Sl  Eldefrawy*  1992;  Farhat  ct  al., 
1994;  Feudd  ct  al.*  2000;  Holden  ct  ah*  1992;  Izhikevich, 
2000),  For  example*  dynamics  of  interspike  time  interval  of 
biological  thermally  sensitive  neurons  with  increasing 
temperature  (which  is  both  its  input  and  bifurcation 
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parameter)*  undergoes  transition  to  chaos  via  period- 
doubling  cascade*  intermittency  and  crises  of  chaotic 
attractors,  emerging  windows  of  periodic  activity*  etc,  as 
shown  in  Fig.  L  (Fcudcl  et  al.*  2000). 

More  typically*  chaotic  neural  dynamics  emerges  at 
macro-level*  in  the  network  of  dynamical  units.  Seminal 
results  of  Freeman  and  co-workers  (Freeman*  1988)  suggest 
that  the  state  of  the  olfactory  bulb  in  olfactory  system*  when 
unperturbed*  is  wandering  within  high-dimensional  chaotic 
attractor.  Applied  input  (odor)  shifts  the  system  to  one  of  its 
low- dimensional  attractors*  'wings’*  that  correspond  to  the 
recognized  odor. 

Dynamics  of  the  network  of  parametrically  coupled 
logistic  maps  was  explored  in  (Farhai*  1997*  2000;  Farhai 
et  al,*  1998),  It  was  shown  that  such  networks  may  have 
enormous  memory  capacity  due  to  the  astronomical  number 
of  different  coexisting  dynamic  attractors.  The  lattices  of 
chaotic  maps  were  studied  in  Dmitriev*  Shirokov*  and 
Starkov  (1997),  Kaneko  and  Tsuda  (2000)  and  Sinha  and 
Ditto  (1998).  Chaotic  dynamics  of  a  network  was  also 
explored  in  Adaebi  and  Aihara  (1997)  and  Hoshino*  Usaba, 
Kashimori*  and  Kambara  (1997),  and  applied  to  an 
engineering  application:  a  search  of  an  optimal  solution  of 
the  traveling  salesman  problem  (Tokuda*  Nagashima*  Sl 
Aihara*  1997).  Encoding  with  the  trajectory  of  the  system's 
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Fig.  I  (a)  Bifurcation  diagram  of  the  modified  Hodgkin-Huxky  model  of  the  thermally  sensitive  neuron.  Interspike  intervals  rrt  versus  bifurcation  parameter, 
temperature  T  (From  Feudei  et  al,,  (2000).  Chaos,  I0<])).  (b)  Bifurcation  diagram  of  quadratic  logistic  map  —  j,),  x,  values  versus  bifurcation 

parameter  R  (c)  Zoomed  in  version  of  (b):  Emerging  of  pcriod-3  orbit, 

state  in  the  phase  space  realized  in  a  network  of  FitzHugh- 
Nagumo  spiking  neurons  was  studied  in  Rabinovich  el  aL 
(2001). 

In  biological  neural  circuits,  input  may  not  be  presented 
by  initial  conditions.  It  is,  rather,  one  of  the  bifurcation 
parameters  (Feudei  el  al.T  2000;  Fukai,  Do*,  Nomura,  & 

Sato,  2000;  Koch,  1999),  Dynamics  of  thermally  sensitive 
neurons  (Braun,  Eckhardt,  Braun.  &  Huber,  2000;  Feudei 
ct  al„  2000;  Gilmore  et  aL,  1999)  is  a  good  illustration  of 
this  idea.  This  input-as-a-bi  furcation -parameter  concept, 
explored  in  Farhat  and  Eldefrawy  (1992),  Farhat  et  al. 

(1994)  and  Holden  et  at.  (1992)  provides  an  insight  of  how 
microscopic  fluctuations  of  an  input  may  be  able  to  change 
the  system's  global  dynamics. 

LL  Why  maps? 

Biologically  realistic  modified  Hodgkin- Huxley  neuron 
models  are  barely  analytically  tractable  due  to  a  huge 
number  of  variables  and  bifurcation  points.  However, 
certain  aspects  of  their  activity,  for  example,  the  dynamics 
of  the  inters  pike  time  intervals,  can  often  be  described  with 
dynamics  of  sine-circle  and  other  maps,  which  are, 
generally,  easier  to  work  with  Ermen  trout  and  Kopell 
(1998)  and  Farhat  and  Eldefrawy  (1992,  1994),  In  this 
article,  as  in  Farhat  (1997,  2000),  Farhat  and  Lee  (1998)  and 
Lee  and  Farhat  (2001)  we  use  the  dynamics  of  quadratic 
logistic  maps  (Figs,  t  and  2,  Eq.  (1))  (Holmgren,  1996;  Ott, 


1993)  as  an  abstract  model  of  a  chaotic  processing  element 

3Cf+j  ~  Rxt(  1  —  Jtj)  (I) 

Bifurcation  diagrams  of  QLM  and  modified  HHM  of 
thermal  neurons  are  shown  in  Fig.  1  (a)  and  (b).  The  reasons 
of  their  striking  resemblance  are  saddle-node,  period¬ 
doubling  and  other  common  basic  bifurcations  which 
underlie  these  dynamics.  Period-doubling  cascade  route  to 
chaos  present  in  both  of  them  is  one  of  the  fundamental 
bifurcation  scenarios  which  is  behind  a  huge  number  of 
processes— from  a  population  dynamics  in  ecological 
systems,  to  chemical  reactions,  like  the  one  of  Belousov- 
Zhabotinsky  (Kaneko  Sc  Tsuda,  2000).  Computational 
abilities  of  the  bifurcation  processes  in  the  logistic  maps' 
dynamics  have  been  studied  extensively  in  a  number  of 
works  (Farhat,  1997,  2000;  Farhat  Sc  Eldefrawy,  1992; 
Farhat  et  al.,  1994;  Lysetskiy  et  al.,  2002),  In  this  article,  we 
focus  on  one  of  the  numerous  bifurcation  processes- 
collapse  of  chaos  to  a  period- 3  orbit  in  the  QLM  dynamics 
and  its  potential  computational  properties. 

1.2.  Etynamics 

Here,  we  briefly  review  the  QLM  dynamics  which  is  used 
in  the  following  sections.  With  bifurcation  parameter  R  < 
3,  the  system  Jt>+)  =/[xf]  (Eq,  (1))  has  a  single  stable  fixed 
point  (Fig,  1(b)),  First  period-doubling  bifurcation  occurs 
when  R  =  3.  With  R  increasing  further,  the  system  under¬ 
goes  period-doubling  cascade  and  at  the  critical  point 
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Fig.  2.  (a)  A  neuron  with  a  sigmoidal  activation  function  as  a  logistic  map  and  its  dynamics,  (b)  Example  of  quadratic  logistic  map  jti+,  = 
dynamics:  period-3  orbit,  (c)  Third  iterate  of  QLM  jp+i  =  =  /[/(fUdll.  R  —  3  75  three  stable  fixed  points  are  about  to  be  bom  via  bifurcations  (d)  The 

distances  A A>  A B>  and  AC  (Eq.  (5)}  to  the  corresponding  bifurcation  points  aw,  br„  and  (Eq.  (3))  are  identical  only  for  a  single  critical  value  of  R  {in  this 
working  range  of  ff). 


Rt  ■*=  3.57  it  becomes  chaotic.  Due  to  the  fractal  structure  of 
the  bifurcation  diagram,  there  is  an  infinite  number  of  values 
of  R,  at  which  chaotic  attractors  the  system  lives  on 
collapse,  producing  stable  periodic  orbits. 

In  this  article,  we  focus  on  the  period-3  orbit  which 
emerges  when  /?3  =  3.828  (Fig.  1(b)  and  (c)).  Its  appear¬ 
ance  is  due  to  three  saddle-node  bifurcations,  giving  birth  to 
three  stable  and  three  unstable  orbits  out  of  chaos.  How  this 
phenomena  happen  can  be  easily  seen  graphically.  Period-3 
orbit  of  the  map/[xf]  (Fig.  2(b))  corresponds  to  a  period- 1 
orbit  (fixed  point)  of  the  map  /3[xr]  =  /(/(/[x,]]J  (Fig.  2(c)). 
Fixed  point  x *  of  the  system  x,+3  =  /3[xf]  can  be  defined  as  a 
point  of  intersection  of  curves  x,+3  —  /3[*,)  and  xJ+3  =  xt 
{Figs.  2(c)  and  3).  Stability  of  x*  is  defined  by  Eq.  (2a-c) 


df\x] 

cLr 


<  i; 


(2a) 


by  Eq.  (3) 

d/3W  ■ 

dx 


(3) 


This  produces  three  neutral  fixed  points  via  three  saddle- 
node  bifurcations.  The  points,  then,  split  into  to  three  stable 
and  three  unstable  solutions.  We  name  them  correspond¬ 
ingly  As,  B Cs  and  #u,  Cu.  The  stable  solutions  exist 
while  they  satisfy  Eq.  (2a),  but  as  R  further  increases  they 
loose  their  stability  via  period- doubling  bifurcations. 


2,  Input-induced  bifurcations 

2  L  Dynamics 


d/3W|_ 

dx  I 


In  order  to  make  the  emergence  of  stable  orbits 
(2b)  compute,  we  shift  function  f*[x]  vertically  with  input  / 
(with  weight  w) 


d/3[*] 

djt 


>  1; 


<2c>  *f+3  “  *1 


(4) 


Fixed  point  x*  is  stable,  neutral  or  unstable  if,  respectively, 
condition  (2a),  (2b)  or  (2c)  is  satisfied. 

When  R  is  slightly  less  than  /?3,  map/3[x]  has  no  stable 
fixed  points  and  its  state  wanders  within  chaotic  attractor 
(Figs.  1(c)  and  2(c)).  However,  if  R  —  R$,  the  curve  x,+3  = 
/3[xJ  touches  the  line  x,+3  =  x,  simultaneously  in  three 
saddle-node  bifurcation  points  asn ,  bcn  and  csn,  defined 


With  this  input,  which  is  an  additional  bifurcation 
parameter,  the  system  demonstrates  quite  different 
dynamics  and  bifurcation  scenario.  With  R  <  it  has 
no  stable  fixed  points  and  lives  on  a  chaotic  attractor. 
Now,  if  we  increase  /  starting  from  /  =  0,  the  curve 
xf+3  -  f[xr]  -  wl  will  touch  the  xl+J  —  xt  line  at  three 
distinct  input  values  (Fig.  3).  This  happens  because  with 
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Fig.  3.  (a)  Function  =  f[xt)  =  /EA/W)J  —  wl  of  QLM3.  R  =  3,805.  (b)  Emergence  and  disappearance  of  fixed  points  via  sequence  of  saddle-node 
bifurcations  with  different  inputs.  Black  and  empty  circles  correspond,  respectively,  to  stable  and  unstable  fixed  poims. 


the  same  R  the  distances 
M=/3[flsn]-fls„  A  B^fi[bsa]-bm 

(5) 

AC=f{ca]-Cn 

at  the  bifurcation  points  am,  btn  and  csn  (defined  by  Eq  (3)) 
are  different  (see  Fig,  2(c)  and  (d)).  The  distances  AA,  A£, 
and  AC  are  defined  as  the  solutions  of  the  equation 

4=/V,*’W  (6) 

with  a  given  R*  and  xm,  Generally,  Eq,  (6)  has  multiple 
solutions.  However,  in  the  range  of  R  we  are  interested  in,  the 
solution  of  Eq.  (6)  with  xm  equal  to  asn,  bsn  and  cm  and  given 
R*  are  unique,  except  for/?  =  ,  when  AA  =  AB  =  AC  =  0 

(Fig,  2(d)), 

Thus,  the  shift  wl  —  A A  induces  a  single  saddle-node 
bifurcation:  the  curve  x,+J  -  wl  touches  the  line 

xr+3  =  xr  at  a  single  point  x,  =  asn.  It  splits  then  into  stable 
fixed  point  A,  and  unstable  Au  (Fig,  3b:  column  A,  /  =  2.0) 
when  the  shift  is  increased.  The  chaotic  attractor  collapses 
and  the  state  of  the  system  converges  to  A,,  as  is  shown  in 
Fig,  4  (/  —  LSI), 

As  /  keeps  increasing,  AE  looses  its  stability  via  period- 
doubling  bifurcation  according  to  Eq,  (2c),  Dynamics  of  x 
undergoes  cascade  of  these  bifurcations  and  becomes 
chaotic  (Fig.  4).  Then,  when  the  shift  wl  =  AiS,  another 


stable  point  #E  emerges  via  saddle-node  bifurcation  (Fig.  3b: 
Column  BJ  —  5,0).  At  this  moment  (chaotic)  dynamics  of  x 
converges  to  the  stable  point  Z?5.  (Figs.  3(b)  and  4),  The 
same  bifurcation  mechanism  underlies  the  emergence  of  the 
stable  fixed  point  Cs  (and  then,  the  loss  of  its  stability)  when 
/  is  negative,  and  the  function  /3[x]  is  shifted  upward. 


<J2‘ 


1  - 


i 


Fig.  4.  Dynamics  of  the  map  jr,+3  - /*[*,!  ”  with  R  *  3.E05,  versus  /. 
Three  input  intervals  are  mapped  onto  three  stable  fixed  poims  A ,  8  and  C. 


M.  Lysetskiy,  JM.  Zurada  /  Neural  Networks  17  (2004)  225-232 


229 


Lei  /In  and  / ^  be  the  inputs  that  produce,  correspondingly, 
a  stable  fixed  point  via  saddle-node  bifurcation  (Eq,  (2a))  and 
loss  of  its  stability  via  period-doubling  bifurcation  (Eq.  (2c)). 
The  map’s  dynamics  at  ihe  critical  points  Ay  B  and  C  can  be 
seen,  then,  as  an  analog  of  a  perceptron,  as  it  divides  the  input 
interval  into  two  subintervals.  One  of  them  produces  a  stable 
fixed  point  and  is  defined  as 

/»  <  M  <  /pd  (7) 

where  sign  of  /  is  the  same  as  of  the  corresponding  distance 
A  (6)+  The  inputs  outside  of  this  interval  do  not  invoke  a 
fixed  point  and  the  system  remains  chaotic  or  stays  on  a 
periodic  orbit. 

It  should  be  noted  that  extension  of  the  system  for  the 
case  whenlandware  vectors  is  rather  straightforward.  The 
QLM3  neuron  xf+3  =/3[rf]  —  maps  multidimen¬ 

sional  space  to  the  induced  stable  fixed  point  in  the  same 
way  as  it  does  with  one-dimensional  input  (Fig.  4). 

2.2.  Simulation 

In  the  following  simulation  example  bifurcation  par¬ 
ameter  is  set  at  R  =  3,808.  The  reason  behind  this  choice  is 
that  the  widest  period-3  orbit  in  the  map's  dynamics 
(Fig.  1(b)  and  (c))  emerges  with  =  3.828;  /?0  =  3.808  is 
slightly  less  than  so  the  system  is  chaotic,  but  its  stable 
states  can  be  produced  by  a  small  input. 

The  saddle-node  bifurcation  points,  at  which  the  curve 
xi+3  =  fl [JCf]  —  w/  touches  Jc,+3  =  x,  line  are  defined  by  Eq, 
(2b):  a*  =0.1604,  <^=0.1501,  6* -0.5157,  k#  = 
0.4851,  —  0.9556,  c ^  =  0,9588.  Input  intervals  Eq.  (7) 

that  produce  emergence  of  fixed  points  A ,  B  and  C  (with 


w  =  0-01)  are  calculated  with  Eq.  (2b  and  c)  as  follows 

A  :  £  =  1.81  <!<!$&  =  3.02 

B:  £  =  4.51  <  /  <  £  =  7.48  (8) 

C :  £  =  -0.53  >  l  >  1%  =  -0.91 

Simulation  results  (Fig.  4)  demonstrate  that  QLM3  neuron 
maps  the  set  of  input  intervals  onto  the  set  of  invoked  stable 
fixed  points  4,  B  and  C. 

It  should  be  noted  that  if  I  &  0  the  map  *f+3  — 
-  wl  no  longer  maps  the  interval  [0,1]  on  itself. 
Two  ‘runaway'  regions  appear  at  the  edges  of  the  interval 
(0,1],  They  are  defined  as;  0  <  x  <  xc  and  1  -  xc  <  x  <  1, 
where  xc  is  the  leftmost  solution  (for  jcc  £  [0,  1])  of 
equation  /3[x*]  -  wl  =  x*  if  the  shift  is  downward 
(/  >  0),  and  xc  is  its  rightmost  solution  (for  xc  £  (0,  1])  if 
/  <  0.  In  our  simulations,  where  the  input  maximum,  /mM  = 
15  is  positive  and  R  =  3.805,  the  critical  value  xc  — 
2.9X  l(T3,  which  slightly  restricts  the  set  of  initial 
conditions  to  [jcc,  I  —  jrJ, 

3,  Learning 

The  sizes  and  centers  of  the  intervals  mapped  by  QLM3 
neuron  to  its  output  states  depend  on  R  (Eqs.  (6)  and  (8),  Fig.  5) 
and  w.  To  adjust  w,  one  of  the  learning  algorithms  for  the 
multilevel  neurons  would  be  appropriate*  for  example,  the 
one  described  in  Malinowski,  Cholewo,  and  Zurada  (1995). 

We  focus  on  the  other  option  provided  by  QLM3 
neuron— the  learning  of  the  bifurcating  parameter  R ,  li 


Fig.  5.  Parameter-depending  (napping:  examples  of  the  input  intervals  inducing  stable  fixed  points  with  different  values  of  bifurcation  parameter  R. 
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Fig.  6.  The  mapping  of  the  input  interval  as  a  function  of  bifurcation  parameter  R. 


is,  in  a  sense,  analogous  to  the  optimizing  the  steepness 
of  the  activation  function  of  a  sigmoidal  neuron  to 
maximize  the  information  of  its  output  (Bell  & 
Sejnowski,  1995). 

The  mapping  of  QLM3  neuron  as  a  function  of  R  is 
shown  in  Fig,  6,  where  the  centers  of  mapped  intervals  (also 
functions  of  R)  are  described  as  curves  GA(R ),  GB(R ), 
GC(R).  To  make  an  input  /0  produce  a  desired  output,  for 
example  a  stable  state  At  R  should  be  adjusted  in  a  way  to 
put  lQ  into  the  interval  mapped  to  A  :  /*(j?)  <  /0  < 

(Eq*  (7)).  Defining  it  more  strictly,  we  want  the  input  to  be  in 
the  center  of  this  interval,  /0  —  GA(R\  defined  with  Eqs.  (5) 
and  (9)  as 

GA(R)  =  j  (£(/?)  +  lUR» 

=  j((f3 [a„,  fi]  -  )  +  (f 3  [a# ,  R]  -  0pd))  (9) 

Thus,  learning  of  mapping  70  to  A  transforms  to  the  task  of 
minimizing  the  distance  %  -  GA(/?)I  (Fig.  6).  The  error 
function  can,  then,  be  defined  as  a  square  of  this  distance 

Eat)  =  do  -  G(R))2  (10) 

Tliis  brings  the  error  gradient  learning  rule  to  the  following 
form 

d  E  . 

SR  =  -c—  =  2 c(/0  -  G(R))(j(R)  (1 1) 

where  C(R)  stands  for  the  derivative  of  G(R)  and  c  is  the 
learning  coefficient. 

To  implement  unsupervised  learning  the  learning 
algorithm  has  to  choose  the  centerline  Gj(R)  which  is 


closest  to  the  applied  input  /0:  E}  -  min,  [£,],;'(  E  A,  R.  C], 
and  then  the  error  Ej  has  to  be  minimized.  In  the  case  of 
supervised  learning  the  corresponding  centerline  is  defined 
by  the  teacher's  signal. 

Examples  of  simulations  of  the  proposed  learning 
algorithm  (Eq,  (11))  are  shown  in  Fig,  7.  The  task  was  to 
leam  mapping  of  input  /0  —  6  to  the  desired  output  B.  With 
initial  conditions  R0  -  3,77  and  RQ  -  3,82  this  input 
produces  chaotic  and  period -2  orbits,  respectively  (Fig.  6). 
After  70  learning  steps  (Eq.  (1 1))  R  converged  to  R  -  3.805 
(Fig,  7(a))  and  it  produced  desired  output  B  (Figs.  4  and  5) 
with  the  errors  E  =  6.7  X  10“4  and  E  -  4.4  X  IG^4  respect- 
ively.  Dynamics  of  E  is  shown  in  Fig.  7(b),  The  learning 
coefficient  used  c  —  2  X  10-6. 

Thus,  this  learning  method  enables  the  QLM3  neuron  to 
adjust  adaptively  its  bifurcation  parameter  R  to  map  certain 
input  intervals  to  specific  stable  states  (classes).  The 
upgrade  of  R  changes  the  mapping  of  all  interval  s/clusters 
(three,  in  our  case)  simultaneously  (Figs,  5  and  6), 


4.  Discussion 

This  article  explores  computational  abilities  of  control¬ 
ling  the  collapse  of  a  chaotic  attractor  to  the  stable  orbits  in 
the  dynamics  of  a  quadratic  logistic  map.  Such  control  is 
implemented  with  an  external  input  (additional  bifurcation 
parameter)  to  the  third  iterate  of  QLM.  The  resulting 
processing  unit,  QLM3  neuron,  demonstrates  a  richer 
repertoire  of  behavior  than  a  classical  artificial  neuron 
with  sigmoidal  activation  function.  Besides  saturated 
regions,  where  inputs  from  certain  interval s/clusters  invoke 
different  stable  slates,  QLM3  neuron  also  produces  chaotic 
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Fig.  7.  (a)  Dynamics  of  the  bifurcation  parameter  R  during  learning,  (b)  Error  E  during  learning. 


or  periodic  dynamics  in  response  to  the  rest  of  the  inputs. 
Another  potentially  useful  property  of  the  QLM3  neuron  is 
its  ability  to  adjust  adaptively  its  mapping  by  learning  the 
bifurcation  parameter  value. 

In  this  article,  we  used  the  emergence  of  period-3  orbit 
out  of  chaos  which  happens  at  R  —  3.828.  However,  the 
appearances  of  a  period-/]  orbit  out  of  chaos,  or,  in  other 
words,  transitions  of  nth  iterate  from  chaos  to  /i  fixed  points 
are  ubiquitous  in  the  quadratic  map's  dynamics.  For 
example,  in  the  range  of  R  from  Rc  «  3.57  (when  the 
map’s  dynamics  first  becomes  chaotic)  to  R  =  4,  there  exist 
(2P  —  2)12 p  windows  of  period  p  orbits,  where  p  is  a  prime 
number  (Ott,  1993). 

The  map  of  /uh  iterate  would  look  structurally  simitar 
to  the  map  of  the  third  iterate  shown  in  Fig.  2(c),  but  it 
would  have  n  saddle-node  bifurcation  points  instead  of  3. 
So,  h  would  be  possible  to  construct  a  QLM  based  neuron 
with  a  large  number  of  stored  stable  states  using  other 
periodic  orbits.  It  might  be  easier,  though,  to  artificially 
generate  a  map  with  required  number  of  saddle-node 
bifurcations. 
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Abstract 

Biologically  inspired  model  of  the  olfactory  cortex  is 
proposed  which  realizes  mapping  of  the  input  pattern 
temporal  structure  to  the  spatial  dynamic  of  the  ensem¬ 
ble  of  output  integrate-and-fire  neurone  The  temporal- 
to- spatial  mapping  and  distributed  representation  of  the 
model  allows  realization  of  both  rough  cluster  classifica¬ 
tion  and  fine  recognition  of  patterns  within  a  duster  in 
parallel  and  as  parts  of  the  same  dynamic  process  The 
temporal  structure  of  the  system  provides  the  frase  for 
the  modeling  of  multi-pattern  processing.  The  model 
is  able  to  extract  components  of  complex  odor  patterns 
(which  are  the  spatio-temporal  sequences  of  neuref  ac¬ 
tivity),  segment  and  bind  them  temporally 


1  Introduction 

Flexible  object  recognition,  feature  binding  and  seg¬ 
mentation,  attention  focusing  and  other  pattern  pro¬ 
cessing  tasks  are  hardly  handled  by  pattern  recogni- 
tion  techniques  based  on  stationary  principles.  On  the 
other  hand,  they  are  successfully  resolved  by  biologi¬ 
cal  neural  systems,  where  different  kinds  of  temporal 
dynamic  and  temporal  correlations  are  believed  to  be 
underlying  principles  [2],[4),(llj. 

Olfaction  is  an  example  of  such  a  system,  in  which 
spatio-temporal  dynamics  was  a  subject  for  numerous 
experimental  studies  [3] ,[9], [10] ,[12], [13],  and  theoreti¬ 
cal  modeling  [l)t  [7],  [16].  However,  most  odor  recogni¬ 
tion  techniques  do  not  make  use  of  temporal  encoding 
and  processing.  In  such  systems  static  patterns  are  rec¬ 
ognized  by  the  stationary  pattern  recognition  methods, 
which  do  not  appear  to  have  much  in  common  with  bi¬ 
ological  temporal  dynamics. 

On  the  other  hand,  there  are  a  number  of  olfactory 


cortex  models  which  are  based  strictly  on  the  known 
biological  data  and  produce  spatio-temporal  dynamics 
similar  to  the  one  of  experimental  data  [1|,  (16).  How¬ 
ever,  they  do  not  explain  the  functional  significance  of 
this  dynamic,  which  is  believed  to  be  related  to  cortical 
information  processing. 

Our  model  is  an  attempt  to  solve  the  problems  of  multi- 
pattern  spatio-temporal  processing  the  brain  is  dealing 
with  using  the  tools  which  are  believed  to  be  present 
in  the  cortex.  The  temporal  structure  of  our  network 
provides  the  base  for  the  solution  of  many  higher-level 
tasks  mentioned  above.  One  of  them  is  the  problem  of 
coarseness-sensitivity  flexibility.  A  course  enough  sys¬ 
tem  cannot  distinguish  fine  variations  of  the  patterns 
within  a  cluster,  on  the  other  hand,  a  sensitive  enough 
system  is  not  able  to  detect  what  cluster  the  slightly 
different  patterns  belong  to.  The  architecture  of  a  net¬ 
work  of  spiking  neurons  described  in  Section  3  is  able 
to  realize  both  fine  and  rough  recognition  of  odors  en¬ 
coded  as  spatio-temporal  sequences. 

Another  group  of  tasks  handled  by  biological  systems 
includes  feature  binding,  segmentation,  attention  fo¬ 
cusing  and  other  multi- pat  tern  recognition  problems. 
There  is  experimental  data  that  suggests  that  in  the 
brain  they  are  solved  with  the  temporal  processing 
[4],  [8],  and  there  are  models  that  propose  the  possi¬ 
ble  mechanisms  [2]  ,[51,(11).  In  section  4  we  show  that 
the  temporal  structure  of  our  model  allows  us  to  realize 
the  multi- pattern  processing  in  the  olfactory  cortex. 


2  Olfactory  System 

An  odor  identity  is  defined  by  a  group  of  physical 
and  chemical  parameters  of  the  odor’s  constituent 
molecules  and  their  relative  concentrations.  However, 
these  parameters  are  not  clearly  determined,  nor  is  the 
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correlation  between  their  candidates  and  the  odor  prop¬ 
erties  [17].  For  the  sake  of  simplicity  we  assume  that 
one  constituent  molecule  possesses  one  of  these  crucial 
parameters  and  corresponds  to  one  of  the  odor  compo¬ 
nents.  So,  the  odors  are  presented  by  the  concentration 
vector  C  —  {cj.ca,  ..Cn},  where  Cj  is  the  concentration 
of  j'-th  molecule. 

In  the  olfactory  systems  odors  are  first  perceived  by 
hundreds  (say  n)  of  receptor  neurons  of  the  olfactory 
epithelium.  These  neurons  are  sensitive  to  different 
kinds  of  molecules  and  respond  selectively  to  their  pres¬ 
ence  with  oscillatory  firing*  An  odor  is  further  encoded 
into  a  periodical  temporal  sequence  of  spatial  patterns 
of  synchronized  oscillatory  neural  activity  of  the  olfac¬ 
tory  bulb/ antennal  lobe  [7],  (10],  [13]*  The  spatial  pat¬ 
terns  are  proved  to  be  correlated  with  the  odor  com¬ 
ponents  [9]  [12],  but  the  functional  significance  of 

their  temporal  structure  is  unclear.  Different  experi¬ 
mental  data  supports  two  major  hypothesis  of  its  pos¬ 
sible  role.  Such  experiments  as  (10),  (13)  suggest  that 
the  temporal  structure  of  firing  of  different  ensembles 
contributes  to  encoding  of  odor  identity  in  a  certain 
combinatorial  way.  Other  data  shows  that  precise  tim¬ 
ing  of  a  spike,  or  the  phase  of  the  periodic  firing,  de¬ 
pends  on  the  stimulus  intensity,  that  is  the  concentra¬ 
tion  of  the  corresponding  component  [3],  [6).  These  two 
hypothesis  could  coexist  and  compliment  one  another 
or  be  the  parts  of  a  more  complicated  neural  code. 

Our  model  realizes  a  mapping  of  the  temporal  relations 
of  input  patterns  into  spatio-temporal  dynamic  of  the 
output  activity.  Although  we  follow  the  Hopfield’s  idea 
[6]  and  assume  that  the  timing  of  a  pattern’s  firing 
encodes  the  concentration  of  an  odor,  our  system  could 
still  be  employed  if  the  temporal  structure  carried  some 
other  functional  significance. 

The  spatio-temporal  patterns  in  the  olfactory  bulb  are 
believed  to  be  formed  in  the  following  way:  the  greater 
the  concentration  of  the  odor  component  applied,  the 
earlier  the  correspondent  ensemble  synchronizes  its  ac¬ 
tivity  and  fires  (7).  According  to  the  Hopfield's  hy¬ 
pothesis  [6J,  the  corresponding  concentrations  Cj  of  n 
constituent  molecules  are  encoded  as  absolute  phases 
of  ensemble's  oscillation,  related  to  ref¬ 
erence  phase  of  the  cortical  oscillation  <j>(rK  The  re¬ 
spective  phase  shifts  are  then  equal: 

The  functional  relation  between  stimulus  intensity  and 
phase  advance  of  the  spikes  has  been  proposed  by  Hop- 
field: 

^  =  aln(^/rf)  (1) 

where  each  phase  shift  is  assumed  to  be  proportional 
to  the  logarithm  of  the  corresponding  concentration,  a 


is  a  coefficient  and  5  is  a  scale  factor  (6). 

Such  logarithmic  scaling  makes  the  relative  phases  of 
spatial  patterns  Invariant  to  different  concentrations  of 
the  same  odor.  The  changing  of  the  concentration  of 
a  multi- component  odor  results  in  a  phase  shift  of  the 
whole  pattern,  while  the  relative  phase-shifts  remain 
constant. 

In  our  model  we  define  an  odor  as  a  cluster  of  odor 
patterns,  that  have  identical  components  with  different 
relative  concentrations.  The  system  has  to  be  rough 
enough  to  he  able  to  recognize  that  different  patterns 
may  belong  to  the  same  duster*  On  the  other  hand,  it 
has  to  be  sensitive  enough  to  distinguish  slightly  differ¬ 
ent  odors  within  a  cluster* 

In  our  network,  recognition  of  an  odor  is  represented 
by  the  firing  of  the  neurons  of  a  specific  ensemble  in  a 
specific  sequence*  Cluster  recognition  and  fine  recog¬ 
nition  are  represented  by  activation  of  different  neural 
sub-ensembles. 


3  The  Model 

Our  network  consists  of  a  layer  of  leaky  integrate- an  d- 
fire  neurons  tij,  interconnected  via  arrays  of  interme¬ 
diate  delay-neurons  dj'  (Figure  1).  The  neurons  uj 
are  also  connected  with  the  layer  of  temporal  inputs 
These  inputs  simulate  activity  of  the  olfactory 
bulb  and  the  neural  layer  functionally  corresponds  to 
the  olfactory  cortex  that  receives  and  processes  those 
patterns. 

The  periodic  inputs  Sj(f),  [j  =  l,„n}  represent  n  spa¬ 
tial  patterns  which  correspond  to  n  components  of  odor 
concentration  vector  C  =  The  inputs  are 

presented  as  follows: 


sj(t)  =  s  ^2  *(t+  <P3  -  kT)  (2) 

fc-i 

where  J(£)  is  Dirac  delta  function,  T  is  the  signal’s 
oscillation  period,  s  is  a  spike’s  amplitude  and  the 
phase  shifts  encode  concentrations  of  constituent 
molecules  according  to  equation  (1),  An  example  of 
input  pattern  is  shown  in  Figure  I* 

There  are  three  types  of  neurons  in  the  layer*  Each  of 
them  is  characterized  by  its  state,  that  is  the  neuron's 
membrane  potential:  Uj(£)  for  the  principal  neurons, 
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for  the  delay  neurons,  and  arj'ft)  for  the  selective 
neurons.  The  neurons  and  inputs  are  connected  with 
weights  w^ur\  uAd<,>  and  (Figure  1) 

As  described  below  in  equation  (3),  a  neuron  u,  re¬ 
ceives  corresponding  input  signal  Sj(t)  from  the  j-th 
input  and  lateral  signals  ij'(t)  (which  are  defined  later 
in  the  section)  from  the  activated  neuron  u* ,  which  are 
propagated  and  delayed  by  the  delay- neurons 


Figure  1:  Network  architecture.  The  neurons  of  acti¬ 
vated  sub-ensemble  fui,dj\uf}  are  shown  in 
bold.  Arrows  and  black  circles  represent,  cor¬ 
respondingly,  excitatory  and  inhibitory  connec¬ 
tions. 

When  a  neuron  receives  a  spike,  its  potential  tij(£)  is 
increased  by  the  weighted  value  if  it  is  a 

spike  from  input  level,  or  w(ll<urJ^l(£)=uitneljr^L[di,(t)J 
(operator  L  will  be  defined  later  m  the  section),  if  it  is  a 
spike  propagated  though  a  delay- neuron  d*k%.  If  the  po¬ 
tential  of  a  neuron  reaches  its  threshold  value  utAre^, 
the  neuron  fires.  Its  output  signal  Ij(l)  =  L|uj{i)|  pro¬ 
duces  a  spike  which  is  propagated  to  the  array  of  delay- 
neurons  that  transfer  the  spike  to  all  other  neurons  in 
the  layer.  At  the  same  time  its  potential  Uj  is  instantly 
reset  to  0  as  shown  in  (4)*  Additionally,  the  poten¬ 
tial  Uj  is  constantly  decreasing  with  decay  coefficient 
k ,  These  mechanisms  are  employed  by  all  the  neurons 
in  the  model. 


The  operator  L  used  above  maps  the  functions  of  a  neu¬ 
ron  membrane  potential  to  the  function  of  the  spikes 
f(t)  which  this  neuron  produces.  So,  for  example, 
^MO]  -  MO  where  /j(£)  is  equal  to  s(t  -  tihr„h) 
and  tlhreth  is  the  time  when  the  value  of  membrane 
potential  u}  reaches  its  threshold. 

The  parameters  of  the  equation  are  set  in  such  a  way 
that  in  order  for  a  neuron  to  fire  it  needs  to  re¬ 
ceive  two  spikes  in  the  narrow  time-window  One 

spike  $j(t),  from  the  corresponding  input  neuron  and 
another,  l{£)  =  from  one  of  the  delay-neurons 

(see  details  in  section  5),  An  exception  Is  made  for  the 
very  first  input  spike  in  the  first  cycle,  which  alone  is 
able  to  activate  the  corresponding  neuron.  This  excep¬ 
tion  is  made  in  order  to  add  to  the  model  the  func¬ 
tional  property  of  the  networks  like  LEGION,  where 
the  global  inhibition  of  the  neurons  depends  on  the 
number  of  the  activated  neurons  (2).  Such  inhibi¬ 
tion  ensures  that  the  fewer  neurons  are  activated,  the 
greater  is  the  probability  for  a  neuron  to  fire.  In  our 
model  the  neurons  are  made  more  sensitive  to  the  very 
first  input  spike,  because  there  is  no  activation  yet  in 
the  neuron  layer.  So,  this  first  spike  is  enough  for  them 
to  fire. 

Delay- neurons  in  the  arrays  are  iutegrate-and-fire  neu¬ 
rons  with  added  inherent  propagation  delays  £>*,  [k  = 
L.m}  defined  as: 


+  £]  wWL[x{<(t)]  (5) 

The  parameters  of  the  equation  make  the  delay- neuron 
work  as  a  delay-operator  (details  in  section  5).  When 
dJk  receives  a  spike  at  time  t  it  fires  at  t  +  Dk*  The 
values  of  delays  Dt t  change  gradually  across  the  array 
as  follows: 


dt 


+ 


-fcuJ(0  +  ^(n<ur)Sj(t)  + 

£  f>(n<ur)*KW] 

*=1 


(3) 


Dk=T{k  ml/2\{k=l,..m}  (6) 


where  T  is  the  oscillation  period  (2),  and  m  is  the  num¬ 
ber  of  the  del  ay- neurons  in  the  array. 


**i(t  )  =  “lhr«h  =*  =  0  (4) 
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As  an  example  we  consider  a  simple  case,  where  an 
odor  with  two  components  {0,  ..tCf ,  cj,  ..0}  is  applied 
with  c/  >  cj.  For  convenience  we  will  be  presenting 


it  further  as  {cr,cj}.  When  this  odor  is  applied,  the 
neuron  uj  receives  the  input  spike  a/(t)  first,  at  the 
moment  t\  and  then  uj  receives  sj{t)  at 

The  single  spike  sj  is  enough  for  uj  to  fire  because  of 
the  exception  mentioned  above.  The  neuron  U[  fires  at 
the  moment  t]  and  sends  spikes  Z,[u/{£)]  to  all  other 
neurons  ,  {j  =  l..n,  j  ^  /}  via  delay-arrays 
The  neurons  receive  the  delayed  signals  from  at 
times  t  ■+  {fc  =  However,  it  is  not  enough 

for  them  to  fire  because  a  spike  from  input  level  is  also 
needed.  Although  all  of  them  show  subthreshold  acti¬ 
vation,  only  the  neuron  uj  which  will  receive  the  input 
spike  sj  will  actually  reach  the  threshold  and  fire.  Fi¬ 
nally  the  neurons  that  fire  are  uTl  uj  and  all  interme¬ 
diate  neurons  d{!  in  the  array  which  connects  them. 


mediate  delay-neuron  indicates  the  relative  con¬ 
centration  of  two  components,  logarithm  of  which  lays 
in  the  vicinity  of  the  delay  DK  of  the  corresponding 
del  ay- n  euron: 


DK-—<a\n 


2L 

Cj 


<dk  + 


X 

2m 


(9) 


The  corresponding  spatio-temporal  dynamic  is  shown 
at  the  Figure  2b,  left  column.  The  parameters  used  in 
the  simulation  are  described  in  Section  5. 


4  Temporal  Segmentation 


The  values  of  Dk  set  by  (6)  ensure  that  one  and  only 
one  of  the  del  ay- neurons  fires  and  sends  a  spike'  to  the 
neuron  uj  within  the  time  window  At^K  Although 
all  of  the  del  ay- neurons  in  the  array  fired,  only  one 
of  them  actually  contributes  to  the  firing  of  the  neu¬ 
ron  uJr  To  distinguish  this  contributing  delay-neuron 
from  others  an  additional  layer  of  neurons  is  added 
(Figure  1),  These  neurons  provide  negative  feedback 

^  to  the  intermediate  neurons  that 

did  not  contribute  to  the  firing  of  uj  (5).  This  is  done 
with  integrate- and- fire  neurons  which  work  as  coinci¬ 
dent  time  neurons  with  the  allowed  time  window  At. 
Their  parameters  and  functional  properties  are  analo¬ 
gous  to  the  neurons  u^. 


~P  =  -*ii(0  +  tu(ne“p>L[d^t| 
+  w(n"‘r)L[u><0] 


So,  which  contributed  to  the  firing  of  uj  stays  un¬ 
changed,  while  the  rest  of  the  intermediate  neurons  are 
suppressed  and  will  not  be  sensitive  to  the  spikes  from 
u*  during  time  of  suppression  ,  which  is  defined  as 
follows; 


m  1  ,  (  w(*up)  \ 

(8) 

Now  at  the  output  level  there  js  a  sequence  of  firing  of 
neurons  u/,  dJK!  and  u/  one  after  another.  Firing  of  the 
neurons  ui  and  uj  indicates  that  an  odor  of  the  duster 
is  recognized.  The  firing  of  the  specific  inter- 


According  to  the  mechanism  of  formation  of  the  input 
temporal  sequence  [7],  if  a  mixture  of  several  odors  is 
applied,  the  corresponding  spatio-temporal  sequences 
are  superimposed  and  the  resulting  input  sequence  con¬ 
tains  patterns  of  all  components  of  each  odor.  The 
stronger  the  component,  the  earlier  its  pattern  fires, 
no  matter  which  odor  the  component  belongs  to. 


□son*’ 

□  ECO 

□  □□□ 

□  EDO 

□  HDD 


fa) 


i  >8  *  kT  i-  l9+«  k+Pf)T 


(b) 


Figure  2:  Temporal  pattern  segmentation,  (a)  Distribu¬ 
tion  of  the  neurons  in  the  layer;  (b)  Spatio- 
temporal  neural  dynamics  during  first  pf  (left 
column)  and  second  pf  {right  column)  cycles. 
Axis  z  corresponds  to  the  membrane  potentials 
ii(f)  or  d(i),  =  1,  ,.p/},  T  is  the  period  of  os¬ 

cillations 


In  order  to  segregate  the  odor  patterns  in  time,  one 
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of  the  ensembles  should  win  and  suppress  the  others 
for  a  period  of  several  cycles.  After  thatt  due  to  the 
neural  fatigue,  the  winning  ensemble  stops  firing  and 
the  second  strongest  ensemble  wins  and  fires  during 
the  next  several  cycles  [2],|5]1(7]t[ll],  To  realize  such 
pattern  segmentation  the  modified  network  from  the 
section  2  with  additional  neural  interaction  and  neural 
fatigue  function  F(p)  is  used: 


+  EX><"'“,>iK,co]  + 

+  E  ^n,<r,ii«.(t)]}F(p)  (io) 

where  a  neuron  ensemble  Elt  [z  -  1  (rrZ}  is  a  group 
of  neurons  that  correspond  to  the  components  of  the 
same  odor-  All  neurons  which  are  not  members  of 
the  same  ensemble  are  interconnected  with  negative 
weights  When  a  neuron  Ui  fires,  it  sends  in* 

hibitory  signals  to  the  neurons  that  do 

not  belong  to  any  of  the  ensembles  Et  in  which  the  neu¬ 
ron  Ui  participates.  The  values  of  the  membrane  poten¬ 
tial  of  those  neurons  are  decreased  by 
They  stay  suppressed  for  the  refractory  period  Tr  dur¬ 
ing  which  they  cannot  fire  regardless  of  the  input  sig¬ 
nals  received*  Tr  is  determined  as  follows: 


_  1  wm*r)s 
”  *  n  -2 wln***w 


(11) 


with  100%  probability.  So,  the  ensemble  which  con¬ 
tains  the  winning  neuron,  or,  in  other  words,  the  odor 
with  the  strongest  component  always  wins  the  compe¬ 
tition  first. 

As  an  example  we  consider  the  case  where  two  odors 
(c/.cj)  and  {cp,  cq}  are  applied  to  the  network  with 
the  following  order  of  the  concentrations  of  their  com¬ 
ponents:  cj  >  cp  >  cj  >  cq.  The  temporal  sequence 
of  input  spikes  is  the  same:{s/(t),  sF(t),  sQ{t)} 
The  first  of  them  sj  (t)  activates  neuron  uj  t  which  sends 
inhibitory  signals  to  up  and  uq,  and  excitatory  signal 
to  uj  via  array  of  the  delay- neurons.  When  input  spike 
sp(t)  appears,  the  corresponding  output  neuron  up  is 
suppressed  and  will  not  respond  to  the  impulse.  Then 
input  sj  activates  neuron  uJt  which  already  received 
excitatory  signal  from  u/.  After  that  sQ(t)  fails  to  ac¬ 
tivate  neuron  uq.  Finally  the  neurons  u/  and  u j  fire, 
while  up  and  uQ  remain  silent.  This  way  the  odor 
{u/,  uj}  is  segmented  from  the  background  and  atten¬ 
tion  m  focused  on  it  for  the  period  of  pf  cycles.  Then, 
due  to  the  neural  fatigue  F[p],  the  neurons  uj  and  uj 
stop  firing  and  up  and  uq  become  activate d,  so  the 
attention  is  now  refocused  at  this  odor*  So,  the  odor 
patterns  are  temporally  segregated  and  processed  one 
at  a  time,  as  is  shown  at  Figure  2. 


5  Simulation 

The  neurons  described  by  (3)  and  (7)  work  as  coin¬ 
cident  time  detectors.  They  fire  if  two  spikes  arrive 
to  a  neuron  within  time  window  At,  the  size  of  which 
is  defined  by  the  parameters  of  the  integrate- and-fire 
neurons  as  follows: 


F(p)  is  a  step  function  with  F[p)  =  1,  if  p  <  pf  and 
F(p)  =  0,  if  p  >  pf ,  The  variable  p  is  the  number 
of  times  the  neuron's  potential  reached  its  threshold, 
and  pf  is  the  number  of  firings  after  which  the  neuron 
becomes  insensitive  to  the  inputs  and  stays  silent  for 
another  p/  cycles. 

The  dynamics  of  the  competition  between  ensembles  is 
a  quite  complicated  process  |2J,  [7],  because  for  a  neu¬ 
ron  in  ensemble,  probability  to  be  suppressed  depends 
on  the  statistical  value  of  the  difference  of  received  in¬ 
hibitory  and  excitatory  spikes. 

In  our  model,  according  to  equation  (10)  and  the  except 
tion  for  the  first  spike  in  the  input  sequence  (Section  3), 
the  first  activated  neuron  suppresses  the  others,  which 
are  not  in  its  ensemble  and  does  not  allow  them  to  fire 


k  Wnetir>s  7 

(!££-')  f‘2> 

where  Ariu}  and  AtflJ  are  the  time  windows  of  neurons 
%  and  x^1  respectively.  In  our  simulation  Uthr*,h  - 
^ihrejhi  so  Af^  —  At^x\  thus  we  will  represent  them 
both  as  At. 

As  well  as  this  is  essentially  the  property  of  the  neu¬ 
rons  used  in  the  model,  there  was  no  need  to  implement 
them  by  the  actual  integrate- and- fire  neurons.  The 
neurons  Uj  were  replaced  by  logic  units  u*,  character* 
ized  by  its  state  uj  (t)  that  has  the  following  properties; 
Without  inputs  u*(t)  is  equal  to  0,  If  u*  receives  two 
spikes  with  weights  u;ftleur?  at  the  moments  and  t2 
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with  h  >  ti»  Then: 

IF  |t2  —  tl|  <  At 

THEN  u*{tf)  =  uth^h 

ELSE  u*(t)  =  0 

The  rules  for  units  xJJt  are  analogous  to  the  rules  of  u* . 
State  of  d*£x  is  defined  as;  d^'(£)  —  tfj[l(£)  =  L(bi(t)j, 

The  parameters  used  in  the  simulation  are:  a  —  4,  6  = 

1T  S  =  1,  U(hr-»h  —  =  rfthreAh  =  L  U^atfuf)  = 

0.75,  =  LI,  =  -8.(3963,  w^int]  =  -39.91, 

k  =  0.2197,  pf  =  3,  Ai  =  5,  n  -  4,  m  =  4tT  -  TR  - 
Ts  =  20. 

Dynamic  of  an  example  simulation  is  shown  in  Figure  2, 
where  the  mixture  of  odors  =  {100,50},  and 

{up,uq}  =  {80,3}  is  applied.  During  first  3  cycles 
the  neurons  of  the  first  odor  fire  in  the  following  se¬ 
quence:  {ti/,^\uj},  (Figure  2b,  left  column).  This 
indicates  that  the  odor  {u/,uj}  is  recognized  with  the 
relative  concentrations  of  its  components  defined  by  (9) 
0  <  ln{”^|  <  1.25.  After  3  cycles  {t£p,d^,ug} 
becomes  activated  and  suppresses  in  its  turn  the  first 
ensemble  (Figure  2b,  right  column).  So,  the  ratio  of 
the  concentrations  of  the  components  of  this  odor  is: 
2.5  <  ln[ft]  <  3.75. 


6  Discussion 

Information  in  the  brain  is  represented  by  spatio- 
temporal  dynamics  of  neural  ensembles  in  a  distributed 
manner.  Precise  phase  or  timing  of  a  single  spike, 
or  some  more  sophisticated  fine  structure  of  the  spike 
trains  are  the  candidates  for  the  brain  neural  codes  [4]. 
Different  kinds  of  temporal  correlation  [2],  [5], [11]  and 
temporal- to- spatial  mapping  [15]  also  proved  to  be  em¬ 
ployed  in  the  biological  pattern  processing. 

So,  the  neural  properties,  principles  and  even  the 
temporal- to- spatial  mapping  used  in  our  model  appear 
to  have  their  biological  prototypes.  Although  our  sys¬ 
tem  does  not  clarify  the  neural  dynamics  of  a  real  bi¬ 
ological  cortex,  we  present  an  attempt  to  understand 
the  principles  of  how  the  brain  could  use  its  “hardware” 
and  temporal  processing  for  the  recognition  tasks. 

It  is  shown  that  the  temporal  structure  of  the  proposed 
system  allows  realization  of  the  recognition  flexibility, 
and  provides  the  base  for  the  realization  of  higher- level 
tasks  of  multi-pattern  processing  such  as  pattern  seg¬ 
mentation,  feature  binding,  attention  focusing  and  oth¬ 
ers. 


AcfaieWed0Ttten£;  This  work  luo*  jporuprcd  in  part  by  the  De¬ 
partment  of  the  Navy,  Office  of  Naval  Research,  Grant  NQ00 
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ABSTRACT 

The  analysis  of  dynamical  systems  with  nonintegrable  con¬ 
tinuous-time  dynamics  is  oftentimes  performed  with  the  help 
of  a  Poincare  map.  This  reduces  the  number  of  observed  di¬ 
mensions  by  one  without  a  significant  loss  of  generality  in 
the  results.  Dynamics  produced  by  the  Poincare  map,  how¬ 
ever,  is  still  a  sequence  of  real-valued  quantities,  difficult  to 
represent  in  terms  of  electronic  signals.  Alternatively,  a  bi¬ 
nary  signal  representing  the  time  intervals  between  consec¬ 
utive  piercings  through  the  Poincare  section  can  be  easily 
implemented  in  electronic  circuitry. 

1.  INTRODUCTION 

We  introduce  a  temporal  sequence  to  be  a  signal  composed 
of  Dirac  deltas  separated  by  a  sequence  of  time  intervals. 
The  choice  of  the  Dirac  delta  representation  has  no  mean¬ 
ing  other  than  marking  certain  instances  in  time.  More  for¬ 
mally  speaking,  a  temporal  sequence  is  a  mapping  Z  R 
(integer-to-real),  which  resembles  the  set-theoretical  defini¬ 
tion  of  a  sequence. 

Temporal  sequences  form  a  metric  space.  This  allows 
for  testing  how  close  two  temporal  sequences  are  to  each 
other.  By  defining  a  suitable  metric,  convergence  of  a  dy¬ 
namical  system  Z  (Z  -+  R)  can  be  tested  The  meaning 
of  the  introduced  dynamical  system  description  is  a  tempo¬ 
ral  sequence  undergoing  iterations.  This  formalism  is  suffi¬ 
cient  to  investigate  Cauchy  convergence  as  well  as  synchro¬ 
nization  of  two  temporal  sequences.  Modifying  the  statis¬ 
tics  of  a  temporal  sequence  as  it  is  being  generated  can  be 
considered  a  form  of  modulation.  Information  inscribed  in 
the  sequence  in  this  manner  would  be  detectable  if  the  statis¬ 
tics  of  the  unmodulated  time  sequence  were  known  when 
received. 

The  main  advantage  of  having  the  signal  in  the  form  of 
a  temporal  sequence  is  its  suitability  for  binary  representa¬ 
tion.  A  simple  electronic  circuit  (alternative  designs  have 
been  introduced  [1])  for  hardware  experiments  with  tempo¬ 
ral  sequences  is  presented  in  subsequent  sections. 

This  work  was  sponsored  by  the  Department  of  Navy,  Office  of  Naval 
Research,  Grant  NOOOl  4*01 -1-0630;  The  contents  of  this  information  does 
not  necessarily  reflect  the  position  of  the  U.S.  government. 


Figure  I:  The  LC  tank  circuit  with  negative  conductance 
G  and  amplitude  regulating  current  > /(«).  The  current- 
controlled  current  source  is  a  Schmitt  trigger  circuit 

2.  SIMPLE  TEMPORAL  SEQUENCE  GENERATOR 

The  circuit  includes  a  simple  harmonic  (linear)  oscillator 
with  a  tank  circuit  and  a  negative  conductance,  as  shown 
in  Fig,  1.  The  negative  conductance  G  supplies  energy  to 
cover  for  the  dissipation  from  the  actual  LC  components. 
In  order  to  insure  bounded  oscillations,  an  amplitude  stabi¬ 
lization  technique  is  used  by  injecting  current  i/(i).  Current 
if  switches  between  one  of  two  constant  levels  70  or  — 70. 
Every  such  switching  decreases  the  magnitude  of  oscilla¬ 
tions.  The  dynamic  equations  of  the  oscillator  represent  a 
second  order  system 

=  -G*-i  +  ij(0  (1) 


Whenever  the  injected  current  is  constant  if  =  ±I0,  the  so¬ 
lution  (v(i),  »(*))  is  linear  about  the  fixed  point  (0,  ±/0),  as 
shown  in  Fig.  2.  The  frequency  of  oscillations  and  the  enve¬ 
lope  of  magnitude  expansion  are  determined  by  the  eigen¬ 
values  of  the  Jacobian  [-°£C  _10/C].  These  eigenvalues 
are  the  complex  conjugate  pair  \  ±  i-jfe.  There¬ 
fore  the  frequency  of  the  oscillations  is  fQ  =  .  and 

the  magnitude  follows  the  envelope  e~&1. 

The  fixed  point  is  a  center  of  the  oscillations  spiraling 
outward.  Current  t  /  selects  the  present  location  of  the  fixed 
point  as  either  positive  or  negative  constant  /0.  The  switch¬ 
ing  occurs  when  the  trajectory  gains  enough  magnitude  to 
reach  the  vicinity  of  the  opposite  fixed  point  location.  Since 
only  a  single  variable  t  is  used  to  determine  the  switching 
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Figure  2:  There  is  only  a  single  fixed  point  in  the  phase 
space,  however  its  Location  is  switched  between  ±J0.  The 
Poincare  section  i  =  —  /,sign(v)  is  shown  in  bold. 


h 

1  V 

1 

-h 

h 

Jo 

Figure  3:  A  Schmitt  trigger  injects  current  if  =  ±J0  to 
switch  the  fixed  point  Location  whenever  inductor  current  t 
reaches  value  ± I *. 


threshold,  the  circuit  can  be  implemented  using  a  simple 
Schmitt  trigger.  If  the  present  fixed  point  is  Iq  and  the  trajec¬ 
tory  coordinate  goes  below  the  opposite  fixed  point  -J0,  the 
fixed  point  is  switched.  Likewise,  if  the  present  fixed  point 
is  -J0  and  the  current  coordinate  exceeds  positive  J0t  the 
fixed  point  is  also  switched.  A  constant  J*  has  been  selected 
to  adjust  the  threshold  point  for  the  fixed  point  switching. 
The  described  switching  occurs  at  an  instance  t  o  as  follows: 


V(tf) 


-T0,  if  t(*o)  =  /. 
To,  if  t(io-)  --la 


(3) 


where  t$  —  ±  A),  The  hysteretic  relationship 

between  the  currents  t  /  and  i,  shown  in  Fig.  3,  is  a  source  of 
nonliearity  of  the  circuit  and  provide  its  ability  to  oscillate 
chaotically. 


b. 


Figure  4:  Graphs  of  functions  /  (top)  and  g  (bottom)*  Note 
that  axes  of  the  graph  in  the  bottom  do  not  intersect  in  the 
origin. 


The  set  of  points  in  the  phase  space,  corresponding  to 
the  change  of  state  of  the  Schmitt  trigger,  is  located  on  the 
bold  discontinuous  curve  shown  in  Fig.  2,  This  curve,  de¬ 
scribed  by  equation  i  =  -J^s ign(u),  is  selected  to  be  the 
Poincare  section  of  the  oscillator  dynamics.  Note  that  only 
variable  v  is  a  meaningful  coordinate  of  the  events  hap¬ 
pening  on  the  curve.  The  Poincare  section  determines  a 
discrete-time  dynamical  system.  It  will  be  useful  to  express 
its  equations  using  the  following  state  description: 

*«+i  =  /(»«)  (4) 

=  j(v«)  (5) 

Sequence  {«„}  represents  the  coordinates  of  consecutive 
intersections  of  the  oscillator  trajectory  with  the  Poincare 
section.  Sequence  {71*}  represents  time  intervals  between 
these  intersections.  Figure  4  shows  graphs  of  functions  / 
and  g  resulting  from  the  oscillator  dynamics  example  with 
Jj/Jo  =  L05  and  parameters  GLC  selected  to  obtain  eigen¬ 
values  A  —  2x(0.05  ±  l.OOi).  With  the  help  of  these  func¬ 
tions  and  equations  (4)  and  (5),  the  dynamics  on  the  Poincare 
section  can  be  easily  calculated  without  having  to  integrate 
the  differential  equations  of  the  oscillator 

Sequences  of  vn  and  Tn  are  not  easy  to  obtain  from  the 
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Figure  5:  Graph  of  50000  intervals  between  the  Schmitt 
trigger  transitions.  Intervals  r*+i  are  shown  (vertically)  in 
terms  of  preceding  intervals  r*  (horizontally). 


Figure  6:  Mapping  Tnk+l  =  Final  intervals  Tnk +1 

are  shown  vertically  versus  preceding  final  intervals  Tnk 
(horizontally).  Counting  from  the  left,  dk+x  is  the  number 
of  the  graph  branch.  Dashed  lines  separate  variable  ranges 
occupied  by  the  branches. 


electronic  implementation  of  the  oscillator  Another  observ¬ 
able,  readily  available  through  measurement,  is  the  state  of 
the  Schmitt  trigger  Let  {r*}  be  the  sequence  of  time  inter* 
vals  between  consecutive  transitions  of  the  Schmitt  trigger 
As  illustrated  in  Fig,  2,  every  interval  r*  is  composed  of  sev¬ 
eral  full  rotations  of  the  trajectory  about  the  fixed  point  and 
one  final  rotation  taking  less  then  75%  of  period  to  hit  the 
bold  line  on  the  opposite  side.  With  the  selection  of  exam¬ 
ple  oscillator  parameters,  the  period  is  To  =  1,  which  can 
be  identified  as  the  maximum  value  of  function  g  m  Fig.  4b. 
All  intervals  Tn  less  then  period  T$  determine  the  final  ro¬ 
tation  before  the  Schmitt  trigger  transitions.  The  ordered 
set  {n  :  Tn  <  T0}  defines  a  sequence  {n*}  of  indexes  of 
such  final  intervals  Tn.  Prior  to  the  n**th  final  interval  Tnk, 
there  are  dk  =  nk  -  nk^}  -  1  full  rotations  T0.  This  pro¬ 
vides  a  relationship  between  time  interval  sequences  {r*} 
and  {Tn}: 

Tk  =  dkTo  +  Tnk  (6) 

A  question  naturally  arises  as  to  whether  sequence  {r*} 
could  be  generated  directly  by  a  mapping.  Even  more  specif¬ 
ically,  is  interval  r*+l  uniquely  determined  by  the  preced¬ 
ing  interval  r*?  A  number  of  interval  pairs  (rk,  rk+\ ),  gen* 
erated  using  Equation  (6),  is  displayed  in  Fig.  5.  The  grid 
pattern  is  an  immediate  consequence  of  the  r-interval  be¬ 
ing  composed  of  a  multiple  of  T0  phis  a  final  interval  Tnk . 
In  fact,  the  final  intervals  carry  all  the  information  regard¬ 
ing  the  r-sequence.  As  shown  in  Fig.  6,  final  interval  Tnk 
uniquely  determines  the  successive  final  interval  T**41. 


Moreover,  the  graph  of  this  relationship  is  composed  of  sev¬ 
eral  branches.  For  a  given  Tn*  ,  the  successive  final  interval 
is  located  on  the  dk+i-$t  branch  of  the  graph.  This  way,  fi¬ 
nal  interval  T„k  uniquely  determines  also  the  interval  r*+ 1* 
The  evolution  of  final  intervals  resembles  iterations  of  the 
Bernoulli  shifts,  A  mapping  h  and  quantization  scheme  q 
can  be  developed  by  interpolating  the  graph  in  Fig.  6: 

Tn>+,  =h{Tn,)  (7) 

dk+l  =q{Tn„)  (8) 

With  the  help  of  these  two  functions  the  mapping  of  r- 
intervals  can  be  explicitly  devised  as 

Tk+ 1  =  l(Tk  mod  r0)T0  +  A(t*  mod  Tb).  (9) 

Given  an  initial  interval.  Equation  (9)  allows  for  generating 
the  r-sequence  same  as  the  one  produced  by  the  original 
differential  equation.  However  if  certain  limited  precision 
level  is  assumed,  mapping  h  will  introduce  uncertainty  in 
the  dynamics.  There  are  infinitely  many  branches  contained 
in  the  graph  of  mapping  h.  At  higher  argument  values,  cor¬ 
responding  to  the  right  side  of  the  graph  in  Fig.  6,  even  a 
slight  fluctuation  of  argument,  resulting  from  physical  re¬ 
ality  of  implementation,  will  result  in  indeterminable  value 
of  the  mapping.  Such  behavior  is  typical  for  nonintegrable 
dynamics.  In  order  to  avoid  relying  on  individual  trajec¬ 
tories  of  the  dynamics  being  inspected,  using  its  statistical 
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characterization  may  prove  beneficial  in  practical  applica¬ 
tions.  Let  D  C  (0,  7q)  be  both  the  domain  and  the  range 
of  mapping  k.  Let  hi  be  the  i-th  branch  of  mapping  h,  and 
unambiguously  h  Ui^i  Domain  D  is  divided  into  in¬ 
finitely  many  segments  Di+  Each  segment  Di  is  a  preimage 
of  the  domain  for  a  certain  branch  £>;  =  h~l(D).  Seg¬ 
ments  D]  through  Du>  separated  by  the  dashed  lines  in  the 
range  of  mapping  are  shown  in  Fig.  6 ,  A  fraction  of 
segment  Di  which  is  mapped  onto  segment  D  j  is  a  subseg¬ 
ment  Assume  the  initial  interval  Tn<x  is  selected 

randomly  from  domain  D ,  If  the  dynamics  determined  by 
Equation  (9)  is  ergodic,  conditional  probability  P^  that  the 
successive  interval  falls  into  segment  Dj  given  that  in¬ 
terval  Tn„  is  in  segment  Dt  equals  the  ratio  of  the  segment 
lengths: 


=  iftrUPj)! 

m 


(10) 


Numbers  Fii  arranged  io  a  matrix  P  =  [Pij]  describe  prob¬ 
abilities  of  transitions  from  segment  Di  to  Dj  at  every  it¬ 
erative  step.  The  first  six  rows  and  columns  of  matrix  P 
evaluated  numerically  read 


0 

0 

0 

0 

0.528 

0.471  . .  ; 

0.121 

0.346 

0.237 

0.137 

0.091 

0.065  ... 

0.114 

0.360 

0.226 

0.146 

0.096 

0.056  ... 

0.133 

0.357 

0.232 

0.136 

0.084 

0.056  ... 

0.093 

0.338 

0.259 

0.148 

0.094 

0.065  ... 

0.121 

0.365 

0.218 

0.135 

0.095 

0.063  ... 

*  •  - 

... 

... 

... 

* . . 

(11) 


Note  that  P  is  the  transition  matrix  of  a  Markov  process 
with  states  dk  associated  with  segments  D i.  Matrix  P  is 
probabilistic,  which  means  that  the  row-sums  are  always 
equal  to  L  If  <f>  is  any  initial  probability  distribution  <£{u)  — 
Pr (do  =  n)r  after  k  iterative  steps  the  distribution  evolves 
to  be  4>Pk .  The  iterations  eventually  approach  the  invariant 
probability  distribution  0P*.  If  mapping  h 

is  mixing  [2]  and  expanding  (|fc'(Tn4r)|  >  l)t  the  invariant 
distribution  is  uniquely  determined  by  the  parameters  of 
the  dynamical  system  [3],  independent  of  the  initial  distri¬ 
bution  The  dynamical  system  (9)  generates  observable 
time  intervals  r*  with  symbols  dk  occurring  according  to 
the  distribution  shown  in  Fig.  7.  This  distribution  is  also 
the  left  eigenvector  of  matrix  P  associated  with  eigenvalue 
equal  I. 

In  conclusion,  the  first-order  statistics  of  the  time  inter¬ 
vals  generated  by  a  dynamical  system  with  oscillatory  dy¬ 
namics  can  be  derived  from  its  binary  observable.  The  se¬ 
quence  of  lime  intervals  which  displays  a  grid  pattern  in  de¬ 
layed  coordinates  contains  a  sequence  of  discrete  symbols. 
Such  a  symbol  sequence  may  be  considered  Markovian  and 
the  frequency  of  inter-symbol  transitions  can  be  estimated 
in  order  to  determine  transition  probabilities. 
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Figure  7:  Invariant  probability  distribution  <j>~  for  the  dy¬ 
namical  system  derived  from  time  intervals  rfc.  Probabili¬ 
ties  of  the  first  12  stales  (Pr(cft  =  1)  through  Pr{dt  =  12)) 
shown. 
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Abstract  -  Quadratic  logistic  map  { QLM)  is  pro¬ 
posed  as  a  generalized  Form  of  an  artificial  neuron 
(AN).  Dynamics  of  the  QLM  not  only  exhibits  com¬ 
putational  abilities,  but  also  has  certain  common  Fea¬ 
tures  with  the  one  of  the  modified  Hodgkin-Huxley 
models  of  a  neuron.  The  rest  state  of  the  QLM  neu¬ 
ron  is  wandering  within  a  chaotic  attractor.  Applied 
input  is  an  additional  bifurcation  parameter  of  the 
system.  Input  of  a  certain  range  induces  emergence  of 
corresponding  stable  orbit.  An  arbitrary  large  num¬ 
ber  of  attractors  can  be  stored  in  a  single  QLM  neu¬ 
ron.  We  explore  the  computational  abilities  of  the 
QLM  dynamics  and  argue  that  it  may  reflect  certain 
aspects  of  dynamics  of  biological  neurons, 

I.  INTRODUCTION 

Limitations  of  the  static  nature  of  artificial  neural  net¬ 
works  (ANN)  have  inspired  investigation  in  neuron  mod¬ 
els  with  inherent  dynamics,  such  as  spiking  and  chaotic 
neurons.  In  networks  made  of  such  neurons  recognition 
is  presented  by  convergence  not  only  to  a  fixed  point, 
but  also  to  a  limit  cycle  or  even  chaotic  attractor  [5],  Al¬ 
though  this  is  a  step  toward  modeling  of  the  real  brain 
processes,  when  these  dynamic  attractors  are  used  in  the 
framework  of  static  NN,  they  simply  replace  static  fixed 
points  with  the  dynamic  ones.  As  a  result,  they  do  not 
bring  principally  new  computational  properties. 

Most  artificial  neurons  commonly  used  for  simulations 
are  oversimplified  versions  of  the  Hodgkin-Huxley  model 
(HHM),  which,  in  turn,  is  a  simplified  model  of  the  real 
neural  processes. 

However,  HHM,  as  well  as  its  modified  versions 
(MHHM),  demonstrate  quite  fascinating  dynamics,  with 
a  variety  of  chaotic  phenomena.  For  example,  dynamics 
of  interspike  time  interval  of  biological  thermally  sensi¬ 
tive  neurons  with  increasing  input  (temperature)  which 
is  a  bifurcation  parameter  of  a  system,  undergoes  transi¬ 
tion  to  chaos  via  period-doubling  cascade,  intermittency 
and  crises  of  chaotic  attractors,  emerging  windows  of  pe¬ 
riodic  activity,  etc.,  as  shown  in  Fig.  I.  [3],  Chaotic  dy- 


3000 


T  (°C) 

Fig.  1.  Bifurcation  diagram  of  the  modified  Hodgkin- Huxley 
model  of  thermally  sensitive  neurons.  Interspike  intervals 
r*  versus  bifurcation  parameter  T  (temperature).  (From 
LLFeudel  et  ah,  ChaosT VoL  10:1,  2000.) 

namics  emerges  in  neural  systems  at  macro  level  as  well. 
Seminal  results  of  W.Freeman  and  co-workers  [II]  show 
that  the  state  of  the  olfactory  bulb,  in  olfactory  system, 
when  unperturbed,  is  wandering  within  high-dimensional 
chaotic  attractor.  This  chaotic  wandering  is  seen  as  a 
solution  search.  Applied  input  (odor)  shifts  the  state  of 
the  system  to  the  one  of  its  low- dimensional  wings,  which 
corresponds  to  the  recognized  odor. 

This  concept  of  chaotic  search  was  further  explored  by 
[I], [8],  and  even  applied  to  an  engineering  application:  a 
search  of  an  optimal  solution  of  the  salesman  problem 
[12]. 

Another  idea  is  that  if  chaotic  wandering  were  con¬ 
trolled,  the  patterns  could  be  stored  as  trajectories  of  the 
system’s  state  in  the  phase  space.  That  would  give  an 
enormous  memory  capacity  and  preserve  some  relations 
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(like  similarity,  for  example)  of  stored  patterns.  However 
tempting,  the  concept  has  not  been  yet  realized. 

Sensitivity  to  small  fluctuations  and  initial  conditions 
is  is  another  application  of  chaos  possibly  employed  by 
the  brain  circuits  to  separate,  for  example,  a  tiny  Input 
from  a  strong  background. 

How  exactly  chaos  is  employed  by  the  brain  and  how 
the  computation  principles  of  biological  neural  circuits 
are  different  from  the  ones  of  ANN  is  not  well  under¬ 
stood.  However,  one  of  the  hints  about  this  difference  is 
the  fact  that  in  biological  systems  an  input  may  not  be 
presented  by  initial  conditions.  Instead,  an  input  can  be 
its  bifurcation  parameter.  It  is  exactly  the  case  in  the 
HHM  and  MHHM,  where  external  current  is  a  bifurca¬ 
tion  parameter  of  the  dynamics  of  the  neuron’s  poten¬ 
tial.  Dynamics  of  thermally  sensitive  neurons,  where  the 
temperature  is  a  bifurcation  parameter  for  a  interspike 
interval  dynamics  [2],[3j,[4]  is  a  good  illustration  of  this 
idea.  This  input- as- a- parameter  hypothesis  could  explain 
how  microscopic  fluctuation  of  an  input  is  able  to  change 
global  dynamics  of  a  system,  which  is  often  the  case  for 
biological  sensory  circuits.  Subthreshold  oscillations  may 
be  the  biological  mechanism  to  push  the  neuron  close  to 
the  bifurcation  state,  and  make  it  sensitive  to  the  changes 
of  bifurcation  parameter  (input)  [7], 


Fig.  2.  Bifurcation  diagram  of  quadratic  logistic  map  = 
/hr*  (I  —  xe).  x 1  values  (vertical  axis)  versus  bifurcation 
parameter  R  (horizontal  axis).  (From  R.  A.  Holmgren, 
Springer-Verlag,  New  York,  1996.) 

Another  hint  about  brain’s  possible  computational 
principles  comes  from  the  necessity  of  a  flexible  dynam¬ 
ical  structure  of  attractors  which  represent  features  and 
objects.  ANN  philosophy  basically  is:  one  pattern  -  one 
attractor.  However,  the  brain  have  to  deal  online  with 
multiple  objects  which  are  composite,  constantly  chang¬ 
ing  in  time  and  with  vaguely  defined  features.  In  order 
to  handle  this  task  the  structure  of  attractors  must  not 


only  be  hierarchical  but  dynamical  and  flexible  as  well. 
Chaotic  attractors  undergoing  structural  transformations 
(as,  for  example,  merging  and  segregation),  caused  by 
changing  of  bifurcation  parameter,  could  be  the  mecha¬ 
nism  used  by  the  brain  for  this  purpose. 

In  this  paper  we  use  the  dynamics  of  quadratic  logis¬ 
tic  maps  (See  Fig.  2  and  Eqn.  (I))  [6j,[I0]  as  an  ab¬ 
stract  model  of  a  chaotic  neuron  with  great  computa¬ 
tional  power.  At  a  first  glance  the  QLM  is  very  far  from 
both  ANN  and  their  biological  prototypes,  but,  surpris¬ 
ingly,  it  turns  out  that  QLM  may  actually  be  one  of  the 
the  missing  links  between  them, 

A  single  classical  neuron  with  a  feedback  and  sigmoidal 
activation  function  xt+1  =  /[xl]  can  be  seen  as  a  logis¬ 
tic  map  with  three  fixed  points  and  their  corresponding 
basins  of  attraction  (Fig.  3).  Sigmoidal  form  of  the  ac¬ 
tivation  function  used  in  ANN  not  only  produces  triv¬ 
ial  dynamics  with  limited  computational  abilities,  if  also 
does  not  reflect  subtle  and  sophisticated  computational 
processes  of  a  real  biological  neuron. 


Fig.  3.  A  neuron  with  a  sigmoidal  activation  function  as  a 
logistic  map. 

Although  transformation  function  of  QLM, 

xt+1  =  itr*(I  -i*)  (I) 

is  even  simpler,  its  dynamics  is  tremendously  rich,  and 
thus,  very  promising  from  a  computational  point  of  view. 
With  bifurcation  parameter  R  increasing  starting  from  a 
low  value  (see  Fig.  2),  the  number  of  stable  fixed  points 
is  constantly  doubling  till  the  system  falls  to  chaotic  at¬ 
tractors.  In  their  turn,  chaotic  attractors  undergo  further 
transformations,  such  as  merging  and  segregation  with 
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other  chaotic  attractors,  collapses  which  produce  stable 
periodic  solutions,  etc* 

Bifurcation  diagrams  of  QLM  and  MHHM  of  ther¬ 
mal  neurons  are  shown  in  Figs.  I.  and  2.  The  reason 
of  their  amazing  resemblance  are  sad  die- node,  period- 
doubling  and  other  common  basic  bifurcations  which  un¬ 
derlie  these  dynamics.  Period-doubling  cascade  route  to 
chaos  present  in  both  of  them  is  one  of  the  fundamen¬ 
tal  bifurcation  scenarios  which  is  behind  a  huge  number 
of  dynamic  processes  -  from  a  population  dynamics  in 
ecological  systems,  to  chemical  reactions,  like  the  one  of 
Belousov-Zhabot  insky  [7]  ,[9]*  So,  it  is  not  an  accident 
that  dynamics  of  logistic  maps  may  qualitatively  reflect 
certain  dynamical  processes  of  biological  neurons. 

The  question  now  is:  what  all  this  variety  of  bifur¬ 
cation  processes  could  give  us  in  terms  of  computation? 
This  is  what  we  focus  on  in  this  paper. 


TI.  DYNAMICS 


When  bifurcation  parameter  R  is  small  (R  <  3),  the 
system  x1^1  =  /[x£]  has  a  single  stable  fixed  point  (Fig. 

2) *  Period-doubling  bifurcation  occurs  when  R  —  3. 
With  R  increasing  further,  system  undergoes  period- 
doubling  cascade  and  at  the  critical  point  Rc  ss  3.57  it 
becomes  chaotic.  The  dynamics  of  the  system  is  quite 
complicated  and  is  still  a  subject  of  research.  We  fo¬ 
cus  on  one  of  the  numerous  bifurcation  phenomena:  "At 
certain  values  of  R ,  for  example  at  R$  zz  3.82$,  chaotic 
attractor  the  system  lives  on  collapses,  producing  stable 
period-3  orbit  {Figs.  2  and  4).  This  happens  due  to  three 
saddle-node  bifurcations,  giving  birth  to  3  stable  and  3 
unstable  orbits  out  of  chaos.  This  period-3  orbit  also  un¬ 
dergoes  a  cascade  of  period-double  bifurcations  and  the 
system  falls  to  chaotic  attractors  again.  Due  to  the  frac¬ 
tal  structure  of  the  bifurcation  diagram  such  windows  of 
periodic  orbits  are  ubiquitous  and  can  be  found  at  any  in¬ 
terval  of  the  bifurcation  parameter.  Period- 3  orbit  of  the 
map  /[x*]  corresponds  to  a  period- 1  orbit  (fixed  point) 
of  the  map  /3[x]  —  /[/(/[x*]]].  Fixed  point  x *  of  the 
system  xt+I  =  /3[xe]  is  defined  graphically  as  a  point  of 
intersection  of  curves  x*+1  —  /3[x£]  and  =  xf  (Fig. 

3) .  x*  is  stable  if: 


df\x] 

dx 


<  I 


(2) 


it  is  neutral  if: 


df{x] 

dx 


=  I 


(3) 


Fig.  4.  Bifurcation  diagram  of  QLM.  Emerging  of  period-3 
orbit.  (From  R. A. Holmgren,  Springer-Veriag,  New  York, 
1996.) 


and  unstable  if: 


dx 


>  1 


w 


With  R  <  Rz  map  /3[x]  has  no  stable  solution  and  thus 
wanders  within  chaotic  attractor  (Figs  2  and  4).  When 
R  =  Rz  the  curve  x£+1  =  /3[x*]  touches  the  line  xt+1  = 
xl  simultaneously  in  three  points  A,  B  and  C\  defined  by 
(5): 


df[x] 


(5) 


which  produce  6  intersections.  This  gives  birth  to  3  un¬ 
stable  and  3  stable  solutions.  We  name  them  correspond¬ 
ingly  Ai,  Bi,  Ci  and  As-  82,  C2  (see  Fig,  5).  Solutions 
A2,  B2,  C2  are  stable  while  they  satisfy  Eqn.  (2).  As 
R  further  increases  they  loose  their  stability  via  period- 
doubling  bifurcation.  In  order  to  make  this  phenomena 
compute,  we  will  explore  another  bifurcation  scenario. 
Let  R  —  R&  be  slightly  less  than  bifurcation  value 
when  there  is  no  stable  fixed  points  and  the  system  lives 
on  the  chaotic  attractor.  Now  we  add  input  wlnp  (w 
stands  for  the  connection  weight)  to  the  system,  which 
shifts  the  function  vertically  in  the  following  way: 


xt+1  =  /V]  -  wlnp  (6) 


Now,  if  we  increase  gradually  Jnp  starting  from  Jnp  —  0, 
the  curve  xt+1  =  /3[x*]  —  wlnp  will  not  touch  the 
xt+I  =  xl  line  at  3  points  simultaneously  as  in  the  pre¬ 
vious  case,  but  at  3  distinct  input  values  (Figs  5  and  6). 
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Fig.  5.  (a)  Function  f{x\  =  /[J[/[x‘]]]  of  QLM,  Ft  =  3.808, 
(b)  Emergence  and  disappearance  of  stable  orbits  via  the 
sequence  of  saddle-node  bifurcations.  Black  and  empty 
circles  correspond,  respectively,  to  stable  and  unstable 
fixed  points. 

This  happens  so  because  the  distances 

A  A  =  f[A]  -  A 
AB  =  f[B]  -  B 

AC  =  f\C]  -  C  (7) 

are  different  and  unique  for  different  points  {~A,B  and  C) 
for  all  R  <  i?3*  as  shown  in  Fig.  6,  These  differences  are 
due  to  the  different  speed  of  the  vertical  changes 
which  depends  on  how  far  the  point  is  from  the  origin. 

When  the  shift  wlnp  -  AA  the  curve  xt+1  =  /3[xc]  - 
wlnp  touches  the  line  xt+1  =  x**  so  a  neutral  fixed  point 
A  appears  via  saddle-node  bifurcation  which  splits  then 
into  unstable  fixed  point  Ai  and  stable  A2  (Fig-  5b).  The 
chaotic  attractor  collapses  and  the  state  of  the  system 
converges  to  the  stable  point  ^2,  as  is  shown  in  Fig,  7 

[7  up  =  1.5). 

As  Jnp  keeps  increasing  two  processes  occur.  Fixed 
point  A2  looses  its  stability  due  to  (4).  Also,  when  the 
shift  wlnp  —  A B,  another  stable  point  B2  emerges  via 
saddle-node  bifurcation  by  the  same  mechanism  (Fig. 
5b).  At  some  moment  the  state  of  the  system  jumps 
from  A2  that  is  looses  its  stability  to  B2  that  becomes 
stable  (Fig,  5b),  The  same  bifurcation  mechanism  un¬ 
derlies  appearance  of  a  stable  fixed  point  C2  when  Inp  is 
negative  and  the  system  is  shifted  upward. 

3  different  input  intervals  have  been  mapped  onto  3 
different  attractors.  By  exploring  windows  with  different 
periods,  any  desirable  number  of  intervals  can  be  stored 
to  a  single  neuron. 


Fig.  6.  Function  f3\x]  ~  /[/[/[x*]]]  of  QLM,  R  =  3.75. 

It  is  worthwhile  to  note  that  the  shape  of  the  function 
/3  \x]  has  not  been  specially  designed  for  this  kind  of  com¬ 
putation,  so  only  3  saddle-node  bifurcation  points  can  be 
explored  with  specific  Ft  —  R&.  It  is  easy  to  construct  a 
function  of  specific  shape  with  numerous  sad  die- node  bi¬ 
furcation  points  which  would  make  the  number  of  input 
intervals  stored  practically  unlimited  with  their  sizes  ar¬ 
bitrary  chosen, 

IIL  SIMULATION 

Simulation  of  the  QLM  dynamics  was  performed  with 
parameter  R^  =  3.808.  The  reason  behind  this  choice 
is  that  the  widest  period-3  orbit  (Fig.  4)  emerges  with 
bifurcation  parameter  ^3  —  3,828,  —  3.808  is  slightly 

less  than  #3,  so  the  system  is  chaotic,  but  its  stable  states 
can  be  produced  by  a  small  input.  The  saddle-node  bi¬ 
furcation  points,  where  the  curve  xi+I  =  /3[x*]  —  wlnp 
touches  the  xt+1  =  xl  line  are  defined  by  (3)  as  follows: 
A  =  0,161*  B  =  0.515,  C  =  0,955. 

Stable  fixed  points  A2,  B2  and  C2  loose  their  stabil¬ 
ity  via  double-period  bifurcations  when  they  reach  the 
following  values:  A\  =  0.150,  B2  —  0.484,  C2  —  0.959, 
defined  by  (4). 

Input  intervals  that  produce  emergence  of  attractors 
A,  B  and  C  with  in  —  I  are  defined,  respectively  as 
follows: 

A  :  H'A  =  0.015  <  wlnp  <  AAb2  =  0.026 

B:  A B  =  0.039  <  wlnp  <  A  S26  =  0.070 

C  :  AC  —  -0.004  >  wlnp  >  A C£  =  -0.007 

(8) 

Simulation  results  are  shown  at  Fig.  7, 
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« - - 1 - i - 1 - 1 - 

-2  Q  2  4  6  Jnp 

Fig,  7,  Dynamics  of  the  map  x1*1  =  /3[ar*]  —  wJnp  (with 
R  =  3,808)  versus  applied  input.  Three  input  intervals 
are  mapped  onto  three  attractors 

IV.  DISCUSSION 

We  have  shown  that  quadratic  logistic  map  can  be  con¬ 
sidered  as  a  more  general  neural  model  with  great  com¬ 
putational  abilities,  QEM  demonstrate  dynamics  which, 
at  certain  level,  has  certain  common  features  with  the 
one  of  biological  neurons. 

Bifurcation  transition  that  we  used  (emerging  of 
period-3  window  out  of  chaotic  attractor) ,  is  only  one  of 
numerous  phenomena  of  QEM  dynamics.  Other  chaotic 
processes  of  QEM  could  provide  us  with  more  compu- 
tationaJ  abilities  and  possibly  with  some  cues  to  under¬ 
stand  certain  aspects  of  dynamics  of  biological  neurons, 
Such  process  as,  for  example,  transformations  of  chaotic 
attractors  of  QEM  may  reflect  dynamic  transformations 
of  the  global  state  of  the  biological  neural  system  during 
the  recognition  process,  Appearance  and  then  merger  of 
several  chaotic  attractors  to  a  single  big  chaotic  attractor 
may  correspond  to  a  recognition  of  parts  of  a  composite 
object,  their  binding,  and  then  recognition  of  the  object 
as  a  whole.  Thus  we  argue  that  dynamics  of  quadratic 
logistic  maps  could  be  a  useful  tool  for  exploring  compu¬ 
tational  abilities  of  artificial  and  biological  neurons, 
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Abstract 

According  to  experimental  studies,  interspike  intervals  (151}  axe  responsible  for  odor  information 
encoding  occurring  in  the  olfactory  bulb.  In  this  paper  a  method  is  described  for  generating  IS]  that 
hypothetically  mimics  the  temporal  encoding  of  odor  information.  The  IS  I  is  generated  by  the  dynamical 
system  derived  from  a  Markov  process  synthesized  to  incorporate  the  IS!  as  its  invariant  distribution. 

The  ISI  distribution  is  the  principal  eigenvector  of  the  transition  matrix  of  Markov  process.  An  algorithm 
using  Nelder  *  Mead  to  synthesize  the  transition  matrix  through  a  cost  function  minimization  is  presented. 

1  INTRODUCTION 

A  sense  on  smell  {or  Olfaction)  allows  vertebrates  and  other  organisms  with  olfactory  receptors  to  identify 
the  odorant  molecules,  A  particular  group  of  receptors  is  excited  by  a  specific  group  of  odorants  (Dickson 
et  al.  1998).  The  sensory  olfactory  receptors  response  to  the  odorants  is  then  passed  to  the  olfactory  bulb 
for  detection.  The  function  of  the  olfactory  bulb  is  to  relay  the  input  obtained  from  the  olfactory  receptors 
into  a  set  of  easily  interpretable  signals  to  the  brain.  The  characteristics  of  the  signals  is  assumed  to  be  the 
following  (Lozowski  et  al,  2004); 

L  The  spikes  mark  instances  of  time  when  the  neurons  fire  and  their  shape  is  not  a  significant  factor. 

2,  The  average  value  of  the  signal  is  affected  by  spikes  since  signal  is  a  time  sequence  of  spikes  and  also 
they  fluctuate  between  occurring  more  or  less  frequently 

3,  The  code  division  of  information  conveyed  by  a  single  signal  is  possible  because  the  interspike  intervals 
may  follow  a  distinct  and  a  repetitive  behavior,  in  simple  words  spikes  may  occur  in  a  certain  temporal 
pattern. 

4,  Synchronization  between  the  signal  sources  may  cause  two  or  more  signals  to  exhibit  cross-correlationdf 
the  synchronized  signals  assume  a  certain  spatial  distribution,  a  set  of  such  signals  will  manifest  a 
spatio-temporal  pattern. 

According  to  a  recently  proposed  principle  of  pattern  processing  with  ISI  distribution  which  hypothetically 
mimics  the  olfactory  system  encoding  (Lozowski  et  al.  2004),  the  feature  of  odor  molecules  could  be  encoded 
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with  1SI  distribution  produced  by  a  dynamical  system.  The  dynamic  system  is  derived  from  a  Markov 
process  synthesized  to  incorporate  the  ISI  distribution  which  is  the  principal  eigenvector  of  Markov  process. 

This  paper  discusses  the  means  of  encoding  odor  information  in  the  olfactory  bulb  and  suggests  an  imple¬ 
mentation  for  encoding  by  using  a  biological  neural  circuit  of  spiking  integrate- and- fire  neurons.  This  paper 
also  discusses  an  algorithm  using  Nelder  -  Mead  to  synthesize  the  transition  matrix  through  a  cost  function 
minimization, 

2  SIGNAL  ENCODING  IN  OLFACTORY  BULB 

Olfactory  receptors  located  in  the  olfactory  epithelium  region  are  hard- wired  to  detect  a  specific  odor  com¬ 
ponent.  On  one  end  of  the  receptor  there  is  a  dendrite  which  receives  signals  from  odorants,  on  the  other  end 
is  an  axon,  which  projects  to  the  olfactory  bulb.  Thus  olfactory  epithelium  the  input  stage  to  the  olfactory 
system. 

The  olfactory  receptors  also  called  sensory  neurons  produce  spikes,  whose  frequency  (/)  of  spiking  depends 
on  the  concentration  of  the  odorants  at  olfactory  epithelium.  The  larger  the  intensity  of  the  odorants  at 
olfactory  epithelium  the  more  frequent  the  sensory  neurons  spike  (White  et  al.  1991).  The  concentration 
information  is  temporally  modulated  at  the  glomerular  inputs  of  the  olfactory  bulb,  therefore  the  perception 
of  odor  intensity  must  be  related  to  the  interspike  intervals  {t  =  Iff)  (Lozowski  et  al.  2004),  If  the  odor 
intensity  increases,  the  intervals  of  spiking  decreases  at  a  different  rate  for  each  basic  odorant  due  to  the 
differences  in  their  conversion  gains. 

A  much  more  compressed  way  to  describe  odors  is  through  distributions  of  interspike  interval  probabilities, 
which  is  more  relevant  to  the  signals  presented  to  the  mitral  inputs.  Let  the  interspike  intervals  be  quantized 
into  N  ranges  with  cutoff  rmax.  Maximum  time  difference  between  evoked  and  spontaneous  activity  of  the 
receptors  is  rmax.  A  single  neural  signal  can  represent  an  odor  with  the  interspike  interval  r  probability 
distribution  a=(ai,a3 ajv), which  is  formally  a  vector  of  probabilities  by  Lozowski  et  al.(2004): 

f <  r  <  ^f7max)  if n<N 

an  =  <  (!) 

[/MtWidz  <  r)  if  n  =  N 

The  transformation  of  the  odor  into  spikes  is  a  dynamical  process.  The  odor  information  is  embedded  in 
the  time  realizations  of  signals  and  can  be  retrieved  only  through  observation  of  these  signals  for  a  period 
of  time.  Statistical  analysis  of  the  neural  signals  retrieves  the  probability  distribution  of  the  interspike 
intervals.  A  simple  stochastic  process  can  be  modeled  to  have  the  statistical  properties  representing  a  given 
odor  through  the  probability  distribution. 

Let  a*  according  to  Lozowski  et  al,(2004)  represent  the  probability  distribution  of  odor  at  a  steady  state, 
i.e.  after  all  transient  response  has  vanished,  A  first  order  approximation  of  a  dynamical  system  for  the  odor 
leads  to  Markov  process  with  an  invariant  probability  distribution  of  a*.  Let  A  (N  x  N)  be  the  transition 
matrix  of  the  Markov  process  (Lozowski  et  al.  2002) 

a(fc  +  1)  ^  Aa(k)  (2) 

The  invariant  distribution  is  the  eigenvector  of  transition  matrix  A  associated  with  the  unit  eigenvalue: 
a*  =  Aa* ,  thus  we  call  the  invariant  distribution  as  the  principal  eigenvector  of  the  Markov  process. 

The  realization  of  the  introduced  Markov  process  is  a  sequence  of  interspike  intervals  rk .  The  interval 
range  is  defined  to  be  Tn  =  ( ^rrmax;rmax)  if  n  <  N,  and  T*  =  [rmoXToo)  otherwise,  where  the  interval 
range  index  n  is  defined  in  the  same  manner  as  in  Rospars  et  al.  (2001).  Operator  A  may  be  developed  by 
using  optimization  to  have  a*  as  its  invariant  distribution  of  interspike  intervals  over  a  period  of  time.  The 


2 


elements  of  the  operator  are  denoted  by  Aij,  so  that  A  —  [*4^],  where  i  and  j  are  the  row  column  indices. 
Number  A^  is  the  probability  that  in  the  Markov  process  (Lozowski  and  Noble  2002)  an  interval  from  the 
range  ry  will  follow  the  interval  from  the  range  r}  (Lozowski  et  al.  2004): 

a  _  Pr(n+ 1  g  Tj  and  rk  g  T,) 
lJ  Pr(rk  €  Ts) 

For  a  given  a  starting  with  some  random  Aj ’s,  Nelder-Mead  Simplex  method  of  optimization  can  be 
used  to  find  the  Aiy  *s  corresponding  to  a  minimum  cost  function  value.  As  all  At/s  are  probabilities,  they 
must  be  real  numbers  in  the  range  (0,1)-  Based  on  this  we  define  a  components  of  cost  function  Eq*  (4) 
namely  quadratic  error  Eq.  (5)  and  penalty  function  Eq.  (6) , 

P(A)  =  Fq(A)  +  Fc(A)  (4) 

Operator  A  is  a  probabilistic  matrix  whose  column  vectors  are  normalized  probability  distributions. 
Therefore,  the  columns  of  A  sum  up  to  1.  The  cost  function  component  Fg(A)  calculates  the  error  of 
deviation  for  this  condition  (Lozowski  et  al.  2004): 

jv  /  N  \  2 

(5) 

The  penalty  function  component  of  the  cost  function  FP(A)  calculates  the  sum  of  individual  element 
penalties,  which  defines  the  error  due  to  a  particular  element  swaying  out  of  range  0  to  1.  This  can  be 
represented  by  the  Eq.  (6), 

n  #  if  Aij  K  0 

Fv(A)  =  il  £  (  My  -  1)  if  Aij  >  1  (6) 

J  0  otherwise 

Operator  A  is  a  well  defined  transition  matrix  of  Markov  process  Eq.  (2)  if  F{A)=  0.  The  goal  of  opti¬ 
mization  procedure  is  to  develop  operator  A  with  the  constraint  that  a*  is  its  principal  eigenvector  associated 
with  eigenvalue  1.  To  simplify  the  operator  synthesis,  matrix  A  will  be  assumed  to  be  diagonal! zable  :  A 

=  BAB~[.  The  diagonal  matrix  A  =  diag(A)  is  composed  of  N  eigenvalues  A=(Ait  A2l . , \N)  of  A.  Let 

Ai  =  1.  The  convergence  of  the  dynamical  system  Eq.  (2)  depends  on  the  radius  of  the  remaining  A/s,  for 
i  >  1,  assuming  that  |A*|  <  r  <  1  and  the  radius  is  kept  low  to  improve  the  convergence  rate*  In  this  paper 
the  value  of  r  was  fixed  to  be  equal  to  0.2,  Operator  A  is  diagonal  in  the  basis  constructed  with  the  column 
vectors  of  B.  Since  A|  =  l,  the  first  column  vector  of  B  is  a\  All  other  entries  j  >  1,  are  the  variables 
in  the  optimization  process.  Their  initial  values  are  selected  randomly  from  the  uniform  distribution  in  the 
range  (-1;!).  The  optimization  process  is  explained  in  the  next  section* 

3  SIMPLEX  METHOD  FOR  FUNCTION  MINIMIZATION 

The  Nelder-Mead  Simplex  method  for  function  minimization  (Nelder  and  Mead  1965)  of  multiple  variables 
is  a  classical,  very  powerful  local  descent  algorithm,  which  depends  on  the  comparison  of  function  values  at 
the  vertices  of  a  general  simplex  and  not  on  the  objective  function  derivatives.  The  comparison  of  function 
values  at  vertices  is  followed  by  the  replacement  of  the  vertex  with  the  highest  value  by  a  new  point  generated 
in  the  process  of  optimization.  The  simplex  adapts  itself  to  the  local  landscape,  and  contracts  on  to  the  final 
minimum* 
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Figure  1 :  An  example  lor  transition  matrix  of  a  four  stage  Markov  process  generated  by  simulator  using  simplex  minlmlza- 
tion  method 


As  described  in  the  previous  section  the  first  column  of  the  matrix  B  is  the  principal  eigenvector  a*,  the 
remaining  N(N  -  1)  variables  in  matrix  B  have  to  be  adjusted  to  obtain  optimum  operator  A.  So,  initially 
considering  minimization  of  a  function  of  N(N  -  1)  variables,  thus  a  simplex  of  N(N  -  I)  dimensions 
needs  N(N  —  1)  +  !  points  which  are  formed  by  randomly  selecting  values  (in  the  range  (-1;1))  for  the 
variables  in  matrix  B  and  corresponding  transition  matrices,  cost  functions  are  formed. 

The  N(N  -  1)  variables  in  each  matrix  are  treated  as  vertices  of  the  simplex,  thus  the  vertices  of  simplex 
are  Pq,Pu  Pn(n-i)-  The  cost  function  at  point  P{  is  denoted  as  Yir  the  highest  value  among  the  cost 
functions  is  denoted  as  Yh  (Tii  1 Vi),  and  the  lowest  value  among  the  cost  functions  is  denoted  as 

YtQri=min?J£  15 Y*)-  P  is  defined  as  the  centroid  of  the  points  with  V*  ^  Yh  and  the  distance  from  point 
Pi  to  Pj  is  represented  by  [PxPj]- 

The  simplex  goes  through  a  sequence  of  geometric  transformations  (exp  ansi  on- step  1,  reflection-step  2S 
contraction-step  6  and  multi-contraction-step  7)  at  each  stage  of  comparison  by  replacing  the  highest  by  a 
new  point.  The  optimization  process  is  as  follows  (R.  Cheiouah  et  aL  2003): 

1.  Determine  the  Y\t  YhyP  and  calculate  P *  =  (1  +  a)P  -  aPh  calculate  the  cost  function  at  Pw  as  Y\ 

2.  If  Y{  <  Yt  then  form  p**  and  go  to  step(3)  else  go  to  step(4)  Pmm  =  (1  +  7 )P*  -  7 P,  calculate  the  cost 
function  at  P**  as  Ymm, 

3.  If  Y*m  <  F*  then  replace  Ph  by  P **  and  go  to  step(5). 

4.  If  Ym  >  Vi  for  all  Yi  except  go  to  step(6). 

5.  Replace  Ph  by  P*  and  go  to  step  (8). 

6.  If  Ym  <  Yh  replace  Ph  by  Pmr  Form  P***  =  0Ph  +  {1  —  fl)P  and  calculate  the  cost  function  at  F**  as 
Y*\ 

7.  If  Y  *  *  >  K/i  then  replace  Pi  by  i  -0:N,  else  replace  Ph  by  P*m 

8.  Check  if  VJ  has  reached  satisfactory  value  or  zero,  if  so  the  optimum  value  for  the  matrix  B  is  /V 

The  transition  matrix  A  shown  in  Fig.  1  is  formed  with  the  matrix  B  obtained  from  simplex  method  of 
minimization.  A  simple  shift-map  can  be  constructed  from  the  probabilistic  operator  by  used  in  approxima¬ 
tion  Eq.  (7),  as  described  in  detail  Pingel  et  ah  (2000), 

f(x)  ~  |  x  —  j  +  Amj  |  +  JV  —iji/j  —  AXJ  <  x  <  j  4-  Aij  —  (7) 

V  M  J  Ml  T^l 
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0  i  x  2  3  4 

Figure  2:  Example  of  a  piecewise  linear  shift  map  /(x)  for  the  transition  matrix  shown  in  Fig,  1 


T,ms 


Figure  3:  An  example  of  spike  train  with  specific  1SI  distribution  produced  in  response  to  the  input  pattern 
fc{0.41, 0.28,0,41, 0,31, 0,34} 


Figure  4:  The  pattern  of  normalised  Input  currents  (a)  and  the  resulting  ISI  distribution  (b)  of  the  spike  train. 
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The  resulting  ISI  Oisfribution  and  spike  train 


T,ms 

spit<e  train  genet sted  with  the  input current  pafiem 


Figure  5:  ISI  distribution  and  spike  train  produced  by  the  ISI  simulator 
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Figure  6:  Transition  matrix  with  principal  eigenvector  a*«{  0.117,0.3755,0,4047,0.1025}  generator  generated  by 
the  simufator 
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Figure  7:  Shlftmap  of  the  transitoln  matrix  shown  in  the  Fig.  6 


A  piece- wise  linear  map  as  shown  in  the  Fig.  2  /  :  [0;  IV]  [0;  N],  [0,  JV]  C  R  is  derived  from  probabilities 
included  in  the  operator  P. 

4  INTERSPIKE  INTERVAL  GENERATOR 

An  example  of  a  spike  train  produced  in  response  to  the  input  pattern  /  —  {0.4  ls  0,28, 0.41  ,  0,31, 0.34}  is 
shown  in  Fig.  3.  The  diagram  of  the  normalized  input  pattern  Inorm  =  {0.24;  0.16;  0.23;  0,18;  0.19}  and  the 
resulting  TSI  distribution  of  10  seconds  spike  train  are  shown  in  Fig.  4.  Here  p(ISIn)  is  the  probability  of 
occurrence  of  the  interspike  intervals. 

Simulator  for  the  interspike  interval  generator  and  can  be  visited  at 
lkttp  :  //ti.uo fl.edu/currentvjork/narem/SIMULATOR/spikegen'.  Fig.  5  shows  the  example  of  a  spike 
train  produced  in  response  to  the  input  pattern  I  =  {0.34,0.55,0.27,  0  44,0.40}  and  the  resulting  ISI  distri- 
bution  of  10  second  spike  train. 

Simulator  for  transition  matrix  generator  can  be  visited  at 
'http  \  ( ld.uofl.edu/cwrrentwork/narem/ SI  MU, LATOR/tranmatnx' .  Transition  matrix  generated  by 
the  simulator  and  the  cost  function  for  the  transition  matrix  with  principal  eigenvector 

a*={0.117, 0.3755, 0.4047, 0.1025}  is  shown  in  the  Fig.  6.  The  resulting  shift  map  of  the  transition  matrix  is 
shown  in  the  Fig.  7. 

5  CONCLUSIONS 

This  paper  describes  a  simplest  method  of  encoding  information  in  temporal  sequences  and  show  input- 
output  interactions  which  lead  to  an  odor  detection  and  encoding  mechanism.  This  paper  also  describes  a 
simulator  which  generates  spike  trains  with  the  input  controlled  distributions  of  interspike  intervals.  This 
paper  also  describes  Simplex  method  for  function  minimization,  a  multi -variable  function  minimization  tech¬ 
nique. 
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Anatomy  of  vertebrate  olfactory  system 
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3  Formal  description  of  neural  signals  4  Temporal  modulation1 

Shape  of  individual  spike  is  insignificant.  Hie  Dirac  delta  function  is  a  model  of  spike  fired  at 


In  ter  spike  interval  densities  6  Stochastic  representation  of  densities 
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Shift  maps  can  generate  interspike  densities  8  Odors  as  shift  maps 
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1 1  Time  sequences  13  Implications 


17  Passible  application 


